Statistik regularisasi masalah terbalik yang salah untuk mereka. Turchin (bagian 1)

Halo, Habr! Hari ini kami ingin memberi tahu Anda apa yang dilakukan laboratorium metode eksperimen fisika nuklir , bagian dari JetBrains Research .

Di mana JetBrains dan di mana fisika nuklir, Anda bertanya. Kami sepakat atas dasar cinta untuk Kotlin, meskipun dalam posting ini kami tidak akan membicarakannya. Kelompok kami berfokus pada pengembangan analisis data, pemodelan, dan perangkat lunak penulisan untuk para ilmuwan, dan karena itu berfokus pada kerja sama dan berbagi pengetahuan dengan perusahaan IT.

Dalam artikel ini kami ingin berbicara tentang metode regularisasi statistik yang kami mempopulerkan , yang diusulkan oleh V.F Turchin pada 70-an abad XX, dan implementasinya dalam bentuk kode dalam Python dan Julia.

Presentasi akan sangat rinci, sehingga mereka yang jelas tentang masalah terbalik dapat langsung menuju ke contoh dan membaca teori dalam artikel ini .

• Terjadinya masalah: mengapa harus ada yang mengatur sama sekali?
• Deskripsi Teoritis tentang Regulasi Statistik
• Contoh praktis tentang:
  
  • Python
  • Julia
• Kesimpulan

Terjadinya masalah: mengapa harus ada yang mengatur sama sekali?

Jika cukup untuk diabaikan, pengukuran apa pun dalam percobaan dapat dijelaskan sebagai berikut: ada perangkat tertentu yang menangkap spektrum atau sinyal dari suatu proses dan menunjukkan beberapa angka sesuai dengan hasil pengukuran. Tugas kita sebagai peneliti, melihat angka-angka ini dan mengetahui struktur perangkat, adalah untuk memahami apa spektrum atau sinyal yang diukur. Yaitu, di muka apa yang disebut masalah terbalik . Jika Anda membayangkan ini secara matematis, kita mendapatkan persamaan ini (yang, kebetulan, disebut persamaan Fredholm dari jenis pertama ):

$$

$$f(y) = \int \limits_a^b dx K(x,y)\varphi(x)$$

Bahkan, persamaan ini menjelaskan hal berikut: perangkat pengukur kami diwakili di sini oleh fungsi perangkat kerasnya $$$K(x,y)$yang bekerja pada spektrum yang dipelajari atau sinyal input lainnya $$$\varphi$menyebabkan peneliti mengamati sinyal output $$$f(y)$. Tujuan peneliti adalah mengembalikan sinyal$$$\varphi$ oleh terkenal $$$f(y)$dan $$$K(x,y)$. Anda juga dapat merumuskan ungkapan ini dalam bentuk matriks, mengganti fungsi dengan vektor dan matriks:

$$

$$f_m = K_{mn}\varphi_n$$

Tampaknya rekonstruksi sinyal bukan tugas yang sulit, karena persamaan Fredholm dan sistem persamaan linear (bahkan overdetermined) memiliki solusi yang tepat. Jadi mari kita coba. Biarkan sinyal yang diukur digambarkan sebagai jumlah dari dua gausses:

$$

$$\varphi(x) = 2*N(2, 0.16) + N(4, 0.04)$$

Sebagai instrumen, kami mengambil integrator yang paling sederhana - sebuah matriks yang menerjemahkan sinyal kami ke jumlah kumulatif menggunakan fungsi Heaviside:

$$

$$K_{mn} = \theta(x_m-y_n)$$

Jenis sinyal yang diukur dan perangkat kami, serta hasil pengukuran ditampilkan pada grafik.

Penting bahwa setiap pengukuran nyata memiliki kesalahan, jadi kami akan sedikit merusak hasil kami dengan menambahkan noise normal, memberikan kesalahan pengukuran lima persen.

Kami akan mengembalikan sinyal dengan metode kuadrat terkecil:

$$

$$\varphi^{} = (K^TK)^{-1}K^Tf$$

Dan sebagai hasilnya kita dapatkan:

Sebenarnya, pada ini kita bisa menyelesaikan artikel, sekali lagi setelah meyakinkan diri kita sendiri tentang ketidakberdayaan metode matematika idealis dalam menghadapi kenyataan fisik yang keras dan kejam, dan pergi untuk merokok setrika.

Tapi pertama-tama, mari kita cari tahu apa yang menyebabkan kegagalan ini terjadi pada kita? Jelas, intinya adalah kesalahan pengukuran, tetapi apa pengaruhnya? Faktanya adalah bahwa Jacques Hadamard (orang yang sama yang menambahkan tanda hubung ke rumus Cauchy - Hadamard) membagi semua tugas menjadi yang benar dan salah.

Mengingat yang klasik: “Tidak masuk akal untuk mencari solusi, jika ada. Ini tentang bagaimana menangani tugas yang tidak memiliki solusi. Ini adalah pertanyaan yang sangat mendasar ... ”- kita tidak akan berbicara tentang tugas yang benar dan segera mengambil yang salah. Untungnya, kita telah bertemu ini: persamaan Fredholm yang ditulis di atas adalah masalah terbalik yang salah - bahkan dengan fluktuasi yang sangat kecil dalam data input (dan bahkan kesalahan pengukuran kami jauh dari sangat kecil), solusi persamaan yang diperoleh dengan cara analitik yang tepat dapat berbeda secara sewenang-wenang dari yang sebenarnya .

Anda dapat membaca bukti pernyataan ini di bab pertama karya klasik akademisi A.N. Tikhonova "Metode untuk memecahkan masalah yang keliru." Buku ini memiliki tips tentang apa yang harus dilakukan dengan tugas yang salah, tetapi teknik yang diuraikan di sini memiliki sejumlah kelemahan yang diperbaiki dalam metode Turchin. Tetapi pertama-tama, kami menjabarkan prinsip umum bekerja dengan tugas yang salah: apa yang harus dilakukan jika Anda menemukan tugas seperti itu?

Karena tugas itu sendiri tidak dapat menawarkan apa pun kepada kita, kita harus melakukan kejahatan kecil: melengkapi tugas dengan data sehingga menjadi benar, dengan kata lain, masukkan beberapa * tambahan informasi a priori tentang tugas * (proses ini disebut regularisasi tugas). Berbeda dengan metode Tikhonov klasik, berdasarkan pengenalan fungsional regularisasi parametrized, metode Turchin memanggil di sisi lain menggunakan metode Bayesian.

Deskripsi Teoritis tentang Regulasi Statistik

Strategi

Kami merumuskan masalah kami dalam hal statistik matematika: sesuai dengan implementasi yang terkenal $$$f$(yang kami ukur dalam percobaan) kami perlu mengevaluasi nilai parameter $$$\varphi$. Fungsional$$$\hat{S}$ menghitung $$$\varphi$berdasarkan $$$f$kami akan memanggil strategi . Untuk menentukan strategi mana yang lebih optimal, kami memperkenalkan fungsi kerugian kuadratik . Fungsi kerugian sebenarnya dapat berupa, mengapa kita memilih yang kuadratik? Karena setiap fungsi kerugian mendekati minimumnya dapat diperkirakan dengan fungsi kuadratik:

$$

$$L(\varphi,\hat{S}[f]) = ||\hat{\varphi}-\hat{S}[f])||_{L_2},$$

Dimana $$$\hat{\varphi}$- solusi terbaik. Kemudian kerugian untuk strategi yang kita pilih ditentukan oleh fungsi risiko :

$$

$$R_{\hat{S}[f]}(\varphi) \equiv E[L(\varphi,\hat{S}[f])] = \int L(\varphi,\hat{S}[f])P(f|\varphi)df.$$

Sini $$$P(f|\varphi)$menentukan kepadatan probabilitas dari ansambel kami, di mana rata-rata kerugian dilakukan. Ensembel ini dibentuk oleh pengulangan beberapa pengukuran hipotetis.$$$f$ untuk diberikan $$$\varphi$. Lewat sini,$$$P(f|\varphi)$ - ini adalah kepadatan probabilitas yang diketahui oleh kami $$$f$diperoleh dalam percobaan.

Menurut pendekatan Bayesian, diusulkan untuk dipertimbangkan$$$\varphi$sebagai variabel acak dengan kepadatan probabilitas apriori $$$P(\varphi)$mengekspresikan keandalan berbagai solusi untuk masalah kita.$$$P(\varphi)$ditentukan berdasarkan informasi yang ada sebelum percobaan. Maka pilihan strategi optimal didasarkan pada meminimalkan risiko posteriori :

$$

$$r_{\hat{S}}(\varphi) \equiv E_{\varphi}E_{f}[L(\varphi,\hat{S}[f])|\varphi]$$

Dalam hal ini, strategi optimal diketahui:

$$

$$\hat{S}[f] = E[\varphi|f] = \int \varphi P(\varphi|f)d\varphi,$$

di mana kepadatan posterior $$$P(\varphi|f)$ditentukan oleh teorema Bayes:

$$

$$P(\varphi|f)= \frac{P(\varphi)P(f|\varphi)}{\int d\varphi P(\varphi)P(f|\varphi)}$$

Pendekatan ini akan memungkinkan kita untuk menentukan varians (fungsi korelasi) dari solusi yang dihasilkan:

$$

$$D(x_1,x_2) = E[\varphi(x_1) - \hat{S}[f](x_1)][\varphi(x_2) - \hat{S}[f](x_2)]$$

Jadi, kami telah mendapatkan solusi optimal untuk masalah kami dengan memperkenalkan kepadatan a priori $$$P(\varphi)$. Bisakah kita mengatakan sesuatu tentang dunia fungsi itu$$$\varphi(x)$yang diberikan oleh kepadatan apriori?

Jika jawaban untuk pertanyaan ini adalah tidak, maka kita harus menerima semua kemungkinan$$$\varphi(x)$sama-sama mungkin dan kembali ke solusi tidak teratur. Karena itu, kita harus menjawab pertanyaan ini dengan tegas.

Inilah yang dimaksud dengan metode regularisasi statistik dalam - regularisasi solusi dengan memperkenalkan informasi a priori tambahan tentang$$$\varphi(x)$. Jika peneliti sudah memiliki informasi apriori (kepadatan apriori$$$P(\vec{\varphi})$), dia hanya bisa menghitung integral dan mendapatkan jawabannya.

Jika tidak ada informasi seperti itu, paragraf berikutnya menjelaskan tentang informasi minimal yang mungkin dimiliki seorang peneliti dan bagaimana menggunakannya untuk mendapatkan solusi yang teratur.

Informasi priori

Seperti yang ditunjukkan oleh para ilmuwan Inggris, di seluruh dunia mereka suka membedakan. Selain itu, jika ahli matematika akan mengajukan pertanyaan tentang legalitas operasi ini, ahli fisika optimis percaya bahwa hukum alam dijelaskan oleh fungsi "baik", yaitu, lancar.

Dengan kata lain, itu membuatnya lebih halus$$$\varphi(x)$lebih tinggi kepadatan probabilitas apriori. Jadi mari kita coba memperkenalkan probabilitas a priori berdasarkan kelancaran. Untuk melakukan ini, kita ingat bahwa pengenalan informasi apriori adalah suatu kekerasan terhadap dunia, memaksa hukum alam untuk mencari cara yang nyaman bagi kita.

Kekerasan ini harus diminimalkan, dan dengan memperkenalkan apriori kepadatan probabilitas, perlu bahwa informasi Shannon di$$$\varphi(x)$terkandung dalam $$$P(\vec{\varphi})$sangat minim. Memformalkan hal di atas, kami memperoleh bentuk kepadatan apriori berdasarkan kelancaran fungsi. Untuk melakukan ini, kami akan mencari ekstrum informasi bersyarat:

$$

$$I[P(\vec{\varphi})] = \int \ln{P(\vec{\varphi})} P(\vec{\varphi}) d\vec{\varphi} \to min$$

Di bawah kondisi berikut:

1. Kondisi kehalusan $$$\varphi(x)$. Biarkan saja$$$\Omega$Merupakan matriks tertentu yang mencirikan kelancaran fungsi. Maka kita mengharuskan nilai fungsional kelancaran tertentu tercapai:
   

   $$

   $$\int(\vec{\varphi},\Omega\vec{\varphi}) P(\vec{\varphi}) d\vec{\varphi} = \omega$$
   Pembaca yang penuh perhatian harus mengajukan pertanyaan tentang menentukan nilai parameter.
   $$$\omega$. Jawabannya akan diberikan lebih lanjut dalam teks.
   
2. Normalisasi probabilitas per unit: $$$\int P(\vec{\varphi}) d\vec{\varphi} = 1$
   Dalam kondisi ini, fungsi berikut akan memberikan fungsionalitas minimum:
   

   $$

   $$P_{\alpha}(\vec{\varphi}) = \frac{\alpha^{Rg(\Omega)/2}\det\Omega^{1/2}}{(2\pi)^{N/2}} \exp(-\frac{1}{2} (\vec{\varphi},\alpha\Omega\vec{\varphi}))$$
   Parameter $$$\alpha$Terhubung dengan $$$\omega$, tetapi karena kami tidak memiliki informasi tentang nilai-nilai spesifik dari kelancaran fungsional, mencari tahu persis bagaimana itu terhubung tidak ada gunanya. Lalu apa yang harus dilakukan dengan$$$\alpha$, Anda bertanya. Di sini tiga cara diungkapkan kepada Anda:
   
   
   1. Nilai parameter yang cocok $$$\alpha$secara manual dan dengan demikian benar-benar melanjutkan ke regularisasi Tikhonov
   2. Rata-rata (mengintegrasikan) dari semua yang mungkin $$$\alpha$dengan asumsi semua mungkin $$$\alpha$sama-sama mungkin
      
   3. Pilih yang paling mungkin $$$\alpha$oleh kepadatan probabilitas posteriornya $$$P(\alpha|\vec{f})$. Pendekatan ini benar karena memberikan perkiraan integral yang baik jika data eksperimen mengandung cukup banyak informasi$$$\alpha$.

Kasus pertama kurang menarik bagi kami. Dalam kasus kedua, kita harus menghitung integral jelek di sini:

$$

$$\left\langle \varphi_i \right\rangle = \frac{\int d\varphi\, \varphi_i P(f|\varphi) \int\limits d\alpha\,P(\alpha) \alpha^{\frac{Rg(\Omega)}{2}} \exp(-\frac{\alpha}{2} (\vec{\varphi},\Omega\vec{\varphi}))}{\int d\varphi P(f|\varphi) \int\limits d\alpha\,P(\alpha) \alpha^{\frac{Rg(\Omega)}{2}} \exp(-\frac{\alpha}{2} (\vec{\varphi},\Omega\vec{\varphi}))}$$

Untuk kasus ketiga, kita dapat memperoleh nilai integral secara analitis untuk kebisingan Gaussian dalam percobaan (ini akan dipertimbangkan melalui bagian).

Perlu juga dicatat bahwa kita belum pernah menggunakannya di mana pun, itu$$$\Omega$Adalah operator kelancaran. Bahkan, kita dapat menggunakan operator lain (atau kombinasi liniernya) di sini, hanya kelancaran fungsi adalah bentuk paling jelas dari informasi apriori yang dapat kita gunakan.

Contoh

Kami berbicara tentang fungsi, tetapi perangkat nyata tidak dapat mengukur tidak hanya sebuah kontinum, tetapi bahkan satu set poin yang dapat dihitung. Kami selalu melakukan pengukuran dalam himpunan titik yang terbatas, oleh karena itu kami terpaksa melakukan prosedur diskritisasi dan transisi dari persamaan integral ke matriks. Dalam metode regularisasi statistik, kami melanjutkan sebagai berikut: kami akan membusuk$$$\varphi(x)$ lebih dari beberapa sistem fungsi $$$\{T_n\}$:

$$

$$\varphi(x) = \sum \limits_n \varphi_n T_n(x).$$

Dengan demikian, koefisien ekspansi ini membentuk beberapa vektor $$$\vec{\varphi}$yang merupakan vektor dalam ruang fungsi.

Sebagai ruang fungsional, kita dapat mengambil ruang Hilbert, atau, misalnya, ruang polinomial. Selain itu, pilihan dasar dalam ruang-ruang ini hanya dibatasi oleh imajinasi Anda (kami mencoba untuk bekerja dengan deret trigonometri Fourier, poligandra dan kubik splines).

Kemudian elemen-elemen dari matriks$$$K$ dihitung sebagai:

$$

$$K_{mn} = (\hat{K}T_n(x))(y_m),$$

Dimana $$$y_m$- titik di mana pengukuran dilakukan. Elemen Matriks$$$\Omega$ kami akan menghitung dengan rumus:

$$

$$\Omega_{ij} = \int\limits_a^b \left(\frac{d^pT_i(x)}{dx}\right)\left(\frac{d^pT_j(x)}{dx}\right)dx,$$

Dimana $$$a$dan $$$b$- batas-batas interval di mana fungsi didefinisikan $$$\varphi(x)$.

Untuk menghitung ulang kesalahan, gunakan rumus dispersi kombinasi linear dari variabel acak:

$$

$$D[\varphi(x)] = D[\sum \limits_n \varphi_n T_n(x)] = \sum\limits_{i,j} \varphi_i\varphi_j cov(T_i(x), T_j(x)).$$

Harus diingat bahwa dalam beberapa kasus, representasi fungsi menggunakan vektor dimensi terbatas menyebabkan hilangnya sebagian atau perubahan informasi. Faktanya, kita dapat menganggap aljabar sebagai sejenis regularisasi, betapapun lemah dan tidak cukupnya untuk mengubah tugas yang salah menjadi tugas yang benar. Tapi, bagaimanapun, kami sekarang telah beralih dari pencarian$$$\varphi(x)$ untuk pencarian vektor $$$\vec{\varphi}$dan di bagian selanjutnya kita menemukannya.

Kasus Gaussian Noise

Kasus ketika kesalahan dalam percobaan didistribusikan menurut Gauss luar biasa
karena solusi analitis untuk masalah kita dapat diperoleh. Karena informasi dan kesalahan apriori memiliki bentuk Gaussian, produk mereka juga memiliki bentuk Gaussian, dan kemudian integral jelek yang kami tulis di atas mudah diambil. Solusi dan kesalahannya adalah sebagai berikut:

$$

$$\vec{\varphi} = (K^T\Sigma^{-1}K + \alpha^*\Omega)^{-1}K^T\Sigma^{-1^{T}}\vec{f}$$

$$

$$\Sigma_{\vec{\varphi}} = (K^T\Sigma^{-1}K+\alpha^*\Omega)^{-1},$$

Dimana $$$\Sigma$- matriks kovarians dari distribusi Gaussian multidimensi, $$$\alpha^*$- nilai parameter yang paling memungkinkan $$$\alpha$, yang ditentukan dari kondisi maksimum kemungkinan kepadatan a posteriori:

$$

$$P(\alpha|\vec{f}) = C'\alpha^{\frac{Rg(\Omega)}{2}}\sqrt{|(K^T\Sigma^{-1}K+\alpha\Omega)^{-1}|}\exp(\frac{1}{2} \vec{f}^T\Sigma^{-1}K^{T}(K^T\Sigma^{-1}K+\alpha\Omega)^{-1}K^T\Sigma^{-1^{T}}\vec{f})$$

Dan jika saya tidak memiliki kesalahan Gaussian?

Bagian kedua dari artikel ini akan dikhususkan untuk ini, tetapi untuk sekarang mari kita uraikan esensi masalahnya.

$$

$$\left\langle \varphi_i \right\rangle = \frac{\int d\varphi\, \varphi_i P(f|\varphi) \int\limits d\alpha\,P(\alpha) \alpha^{\frac{Rg(\Omega)}{2}} \exp(-\frac{\alpha}{2} (\vec{\varphi},\Omega\vec{\varphi}))}{\int d\varphi P(f|\varphi) \int\limits d\alpha\,P(\alpha) \alpha^{\frac{Rg(\Omega)}{2}} \exp(-\frac{\alpha}{2} (\vec{\varphi},\Omega\vec{\varphi}))}$$

Masalah utama adalah bahwa integral yang mengerikan ini, pertama multidimensi, dan kedua, dalam batas tak terbatas. Selain itu, sangat multidimensi, vektor$$$\vec{\varphi}$ dapat dengan mudah memiliki dimensi $$$m = 30-50$, dan metode grid untuk menghitung integral memiliki kompleksitas tipe $$$O(n^m)$, oleh karena itu, tidak berlaku dalam kasus ini. Saat mengambil integral multi-dimensi, integrasi Monte Carlo berfungsi dengan baik.

Selain itu, karena batas kami tidak terbatas, kami harus menggunakan metode pengambilan sampel yang penting, tetapi kemudian kami harus memilih fungsi untuk pengambilan sampel. Untuk membuat semuanya lebih otomatis, Anda harus menggunakan Markov Chain Monte Carlo (MCMC) , yang dapat secara mandiri mengadaptasi fungsi pengambilan sampel ke integral yang dihitung. Kami akan berbicara tentang penerapan MCMC di artikel selanjutnya.

Bagian praktis

Implementasi pertama dari metode regularisasi statistik ditulis kembali pada tahun 70-an di Algol dan berhasil digunakan untuk perhitungan dalam fisika atmosfer. Terlepas dari kenyataan bahwa kami masih memiliki sumber-sumber algoritma tulisan tangan, kami memutuskan untuk menambahkan sedikit modernisme dan membuat implementasi dengan Python, dan kemudian pada Julia.

Python

Instalasi

Instal melalui pip:

pip install statreg

atau unduh kode sumber .

Contohnya

Sebagai contoh, pertimbangkan cara menggunakan modul stareguntuk memulihkan data untuk persamaan matriks dan integral.

Kami mengimpor paket-paket ilmiah yang diperlukan.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.integrate import quad
%matplotlib inline

Kami menentukan sinyal yang benar, yang akan kami pulihkan.

a = 0
b = 5
#  
phi = lambda x: 4*norm.pdf(x-2, scale=0.4) + 2*norm.pdf(x-4, scale = 0.5)
x = np.linspace(a, b,100)
plt.plot(x, phi(x));

Tentukan kernel dan operasi konvolusi fungsi (Catatan: np.convolutionkhusus untuk array):

kernel = lambda x,y : np.heaviside(x-y, 1) #  
convolution =  np.vectorize(lambda y: quad(lambda x: kernel(x,y)*phi(x), a,b)[0])

Kami menghasilkan data yang diukur dan mengeluarkannya menggunakan distribusi normal:

y = np.linspace(a, b, 50)
ftrue = convolution(y)
sig = 0.05*ftrue +0.01 #  
f = norm.rvs(loc = ftrue, scale=sig)
plt.errorbar(y, f, yerr=sig);

Kami memecahkan persamaan integral

Kami mengimpor kelas solver dan auxiliary untuk diskritisasi:

from statreg.model import GaussErrorUnfolder
from statreg.basis import CubicSplines

Sebagai dasar fungsional untuk diskritisasi, kami menggunakan splines kubik, dan sebagai kondisi tambahan, kami mengindikasikan bahwa fungsi tersebut mengambil nilai nol pada tepinya.

basis = CubicSplines(y, boundary='dirichlet')
model = GaussErrorUnfolder(basis, basis.omega(2))

Kami memecahkan persamaan:

phi_reconstruct = model.solve(kernel, f, sig, y)

Kami sedang menyusun jadwal solusi:

plt.plot(x,phi(x))
phir = phi_reconstruct(x)
phiEr = phi_reconstruct.error(x)
plt.plot(x, phir, 'g')
plt.fill_between(x, phir-phiEr, phir + phiEr, color='g', alpha=0.3);

Kami memecahkan persamaan matriks

Kami mengimpor kelas solver dan auxiliary untuk diskritisasi:

from statreg.model import GaussErrorMatrixUnfolder
from statreg.basis import CubicSplines

Untuk mendapatkan matriks, kami menggunakan basis fungsional kami, tetapi jelas bahwa matriks dapat diperoleh dengan cara apa pun.

cubicSplines = CubicSplines(y, boundary='dirichlet')
omega = cubicSplines.omega(2)
Kmn = cubicSplines.discretizeKernel(kernel,y)

Kami memecahkan persamaan matriks:

model = GaussErrorMatrixUnfolder(omega)
result = model.solve(Kmn, f, sig)

Buat bagan:

phir = lambda x: sum([p*bf(x) for p, bf in zip(result.phi,cubicSplines.basisFun)])
plt.plot(x,phir(x))
plt.plot(x,phi(x));

Julia

Seperti yang kami sebutkan, pengembangan lebih lanjut dari teknik ini membutuhkan integrasi Monte Carlo yang canggih. Kita dapat menggunakan beberapa jenis modul dengan Python (misalnya, kita bekerja dengan PyMC3), tetapi kita, antara lain, berpartisipasi dalam proyek bersama dengan Max Planck Institute di Munich.

Proyek ini disebut Bayesian Analysis Toolkit . Tujuannya adalah untuk menciptakan kerangka kerja dengan alat untuk metode analisis Bayesian, terutama termasuk alat untuk MCMC. Sekarang tim sedang mengerjakan versi kedua kerangka kerja, yang ditulis dalam Julia (yang pertama ditulis dalam C ++ buruk). Salah satu tugas kelompok kami adalah menunjukkan kemampuan kerangka kerja ini dengan menggunakan contoh regularisasi statistik, jadi kami menulis sebuah implementasi di Julia .

using PyCall
include("../src/gauss_error.jl")
include("../src/kernels.jl")

a = 0.
b = 6.

function phi(x::Float64)
    mu1 = 1.
    mu2 = 4.
    n1 = 4.
    n2 = 2.
    sig1 = 0.3
    sig2 = 0.5

    norm(n, mu, sig, x) = n / sqrt(2 * pi*sig^2) * exp(-(x - mu)^2 / (2 * sig^2))
    return norm(n1, mu1, sig1, x) + norm(n2, mu2, sig2, x)
end
x = collect(range(a, stop=b, length=300))

import PyPlot.plot

myplot = plot(x, phi.(x))
savefig("function.png", dpi=1000)

Kali ini kami menggunakan inti yang berbeda, kami tidak akan mengambil langkah integrasi, tetapi konvolusi dengan Gaussian, yang sebenarnya mengarah pada "blur" tertentu ke data kami:

function kernel(x::Float64, y::Float64)
    return getOpticsKernels("gaussian")(x, y)
end

convolution = y -> quadgk(x -> kernel(x,y) * phi(x), a, b, maxevals=10^7)[1]
y = collect(range(a, stop = b, length=50))
ftrue = convolution.(y)
sig = 0.05*abs.(ftrue) +[0.01 for i = 1:Base.length(ftrue)]
using Compat, Random, Distributions
noise = []
for sigma in sig
    n = rand(Normal(0., sigma), 1)[1]
    push!(noise, n)
end
f = ftrue + noise
plot(y, f)

Demikian pula, kami mengambil dasar splines dengan ujung tetap:

basis = CubicSplineBasis(y, "dirichlet")
Kmn = discretize_kernel(basis, kernel, y)
model = GaussErrorMatrixUnfolder([omega(basis, 2)], "EmpiricalBayes", nothing, [1e-5], [1.], [0.5])
result = solve(model, Kmn, f, sig)
phivec = PhiVec(result, basis)

x = collect(range(a, stop=b, length=5000))
plot(x, phi.(x))

phi_reconstructed = phivec.phi_function.(x)
phi_reconstructed_errors = phivec.error_function.(x)

plot(x, phi_reconstructed)
fill_between(x, phi_reconstructed - phi_reconstructed_errors, phi_reconstructed + phi_reconstructed_errors, alpha=0.3)

Contoh dunia nyata

Sebagai contoh analisis data nyata, kami akan mengembalikan spektrum hamburan elektron dari campuran hidrogen-deuterium. Dalam percobaan, spektrum integral diukur (yaitu, jumlah elektron di atas energi tertentu), dan kita perlu mengembalikan spektrum diferensial. Untuk data ini, spektrum awalnya direkonstruksi menggunakan fitting, sehingga kami memiliki dasar untuk memeriksa kebenaran algoritma kami.

Beginilah tampilan spektrum terintegrasi awal:

Dan - hasil restorasi:

Analisis dengan fit memiliki tiga kelemahan utama:

• , .
• , .
• .

Regulasi statistik menghindari semua masalah ini dan memberikan hasil model-independen dengan kesalahan pengukuran. Solusi yang diperoleh dengan regularisasi sesuai dengan kurva pas. Perhatikan dua puncak kecil pada 25 dan 30 eV. Diketahui bahwa puncak pada 25 eV terbentuk selama hamburan ganda, dan dipulihkan oleh pas, karena jelas ditentukan dalam fungsi pemasangan. Puncak 30 eV mungkin anomali statistik (kesalahannya cukup besar pada titik ini), atau, mungkin, menunjukkan adanya hamburan disosiatif tambahan.

Kesimpulan dan pengumuman bagian selanjutnya

Kami memberi tahu Anda tentang teknik yang berguna yang dapat disesuaikan dengan banyak tugas analisis data (termasuk pembelajaran mesin), dan mendapatkan jawaban "jujur" yang jujur ​​- solusi paling rasional untuk persamaan dalam menghadapi ketidakpastian yang disebabkan oleh kesalahan pengukuran. Sebagai bonus yang menyenangkan, kami mendapatkan nilai untuk kesalahan keputusan. Mereka yang ingin berpartisipasi dalam pengembangan atau menerapkan metode regularisasi statistik dapat berkontribusi dalam bentuk kode dengan Python, Julia atau pada hal lain.

Pada bagian selanjutnya kita akan berbicara tentang:

• Menggunakan MCMC
• Dekomposisi Cholesky
• Sebagai contoh praktis, kami mempertimbangkan penggunaan regularisasi untuk memproses sinyal dari model detektor orbital proton dan elektron
  

Referensi

• Karya asli Turchin
• Disertasi tentang penerapan regularisasi dalam fisika atmosfer
• Publikasi kami tentang penerapan regularisasi statistik
• Implementasi python
• Implementasi pada Julia
• Implementasi pada ALGOL-60 , untuk pecinta arkeologi
• Perangkat analisis Bayesian

Diposting oleh Mikhail Zeleny , Peneliti di Laboratorium Metode Percobaan Fisika Nuklir di JetBrains Research .