👆 👦🏽 🏄 Statistische Regularisierung falscher inverser Probleme zu ihnen. Turchin (Teil 1) 👨🏻‍🔧 ⛸️ 👩🏾‍💼

Hallo Habr! Heute möchten wir Ihnen sagen, was das Labor für kernphysikalische Experimentiermethoden , das Teil von JetBrains Research ist , tut .

Wo sind JetBrains und wo ist Kernphysik? Wir waren uns aus Liebe zu Kotlin einig, obwohl wir in diesem Beitrag nicht über ihn sprechen werden. Unsere Gruppe konzentriert sich auf die Entwicklung von Datenanalyse-, Modellierungs- und Schreibsoftware für Wissenschaftler und konzentriert sich daher auf die Zusammenarbeit und den Wissensaustausch mit IT-Unternehmen.

In diesem Artikel möchten wir über die von uns populär gemachte Methode der statistischen Regularisierung sprechen , die von V. F. Turchin in den 70er Jahren des 20. Jahrhunderts vorgeschlagen wurde, und über ihre Implementierung in Form von Code in Python und Julia.

Die Präsentation wird sehr detailliert sein, sodass diejenigen, die sich über die inversen Probleme im Klaren sind, direkt zu den Beispielen gehen und die Theorie in diesem Artikel lesen können .

Auftreten des Problems: Warum sollte jemand überhaupt regulieren?

Wenn es ausreicht, zu ignorieren, kann jede Messung im Experiment wie folgt beschrieben werden: Es gibt ein bestimmtes Gerät, das das Spektrum oder Signal eines Prozesses erfasst und einige Zahlen gemäß den Messergebnissen anzeigt. Unsere Aufgabe als Forscher, diese Zahlen zu betrachten und die Gerätestruktur zu kennen, ist es, das gemessene Spektrum oder Signal zu verstehen. Das heißt, angesichts des sogenannten inversen Problems . Wenn Sie sich das mathematisch vorstellen, erhalten wir diese Gleichung (die übrigens die Fredholm-Gleichung der ersten Art genannt wird ):

f (y) = \int_{a}^{b} d x K (x, y) φ (x)

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

Tatsächlich beschreibt diese Gleichung Folgendes: Unser Messgerät wird hier durch seine Hardwarefunktion dargestellt

K (x, y)

$K(x,y)$ , das auf das untersuchte Spektrum oder ein anderes Eingangssignal einwirkt

φ

$\varphi$ , wobei der Forscher das Ausgangssignal beobachtet

f (y)

$f(y)$ . Ziel des Forschers ist es, das Signal wiederherzustellen

φ

$\varphi$ bekannt

f (y)

$f(y)$ und

K (x, y)

$K(x,y)$ . Sie können diesen Ausdruck auch in Matrixform formulieren und Funktionen durch Vektoren und Matrizen ersetzen:

f_{m} = K_{m n} φ_{n}

$f_m = K_{mn}\varphi_n$

Es scheint, dass die Signalrekonstruktion keine schwierige Aufgabe ist, da sowohl die Fredholm-Gleichung als auch das (sogar überbestimmte) lineare Gleichungssystem eine genaue Lösung haben. Also lass es uns versuchen. Das gemessene Signal sei als die Summe zweier Gausse beschrieben:

φ (x) = 2 * N (2, 0.16) + N (4, 0.04)

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

Als Instrument nehmen wir den einfachsten Integrator - eine Matrix, die unser Signal mithilfe der Heaviside-Funktion in eine kumulative Summe umwandelt:

K_{m n} = θ (x_{m} - y_{n})

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

Die Art des gemessenen Signals und unseres Geräts sowie das Messergebnis sind in der Grafik dargestellt.

Es ist wichtig, dass jede echte Messung einen Fehler aufweist, daher werden wir unser Ergebnis ein wenig verderben, indem wir normales Rauschen hinzufügen, was einen Messfehler von fünf Prozent ergibt.

Wir werden das Signal nach der Methode der kleinsten Quadrate wiederherstellen:

φ^{М Н К} = (K^{T} K)^{- 1} K^{T} f

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

Und als Ergebnis bekommen wir:

Eigentlich könnten wir damit den Artikel beenden, nachdem wir uns erneut von der Hilflosigkeit idealistischer Methoden der Mathematik angesichts der harten und rücksichtslosen physikalischen Realität überzeugt haben, und Lötkolben rauchen.

Aber lassen Sie uns zuerst herausfinden, warum uns dieses Versagen passiert ist. Natürlich geht es um Messfehler, aber wie wirken sie sich aus? Tatsache ist, dass Jacques Hadamard (derselbe, der der Cauchy-Hadamard-Formel einen Strich hinzugefügt hat) alle Aufgaben in richtig gestellte und inkorrekte Aufgaben unterteilt hat .

Erinnerung an die Klassiker: „Es macht keinen Sinn, nach einer Lösung zu suchen, wenn überhaupt. Es geht darum, wie man mit einer Aufgabe umgeht, für die es keine Lösung gibt. Dies ist eine zutiefst grundlegende Frage ... “- Wir werden nicht über die richtigen Aufgaben sprechen und sofort die falschen aufgreifen. Glücklicherweise haben wir dies bereits erreicht: Die oben geschriebene Fredholm-Gleichung ist ein falsches inverses Problem - selbst bei unendlich kleinen Schwankungen der Eingabedaten (und selbst unsere Messfehler sind alles andere als infinitesimal) kann die Lösung der Gleichung, die auf exakte analytische Weise erhalten wird, willkürlich von der wahren abweichen .

Sie können den Beweis dieser Aussage im ersten Kapitel der klassischen Arbeit des Akademikers A.N. Tikhonova "Methoden zur Lösung schlecht gestellter Probleme." Dieses Buch enthält Tipps, wie Sie mit falschen Aufgaben umgehen können. Die dort beschriebene Technik weist jedoch eine Reihe von Nachteilen auf, die bei der Turchin-Methode behoben wurden. Aber zuerst beschreiben wir die allgemeinen Prinzipien der Arbeit mit falschen Aufgaben: Was tun, wenn Sie auf eine solche Aufgabe stoßen?

Da die Aufgabe selbst uns nichts bieten kann, müssen wir ein kleines Verbrechen begehen: Ergänzen Sie die Aufgabe mit Daten, damit sie korrekt werden. Geben Sie also einige * zusätzliche a priori Informationen über die Aufgabe * ein (dieser Vorgang wird als Regularisierung der Aufgabe bezeichnet). Im Gegensatz zur klassischen Tikhonov-Methode, die auf der Einführung parametrisierter Regularisierungsfunktionen basiert, werden bei der Turchin-Methode dagegen Bayes'sche Methoden verwendet.

Theoretische Beschreibung der statistischen Regularisierung

Strategie

Wir formulieren unser Problem in Form einer mathematischen Statistik: nach einer bekannten Implementierung

f

$f$ (die wir im Experiment messen) müssen wir den Wert des Parameters bewerten

φ

$\varphi$ . Funktionell

\hat{S}

$\hat{S}$ rechnen

φ

$\varphi$ aufgrund

f

$f$ Wir werden Strategie nennen . Um festzustellen, welche Strategien optimaler sind, führen wir eine quadratische Verlustfunktion ein . Die reale Verlustfunktion kann eine beliebige sein. Warum wählen wir die quadratische? Weil jede Verlustfunktion nahe ihrem Minimum durch eine quadratische Funktion angenähert werden kann:

L (φ, \hat{S} [f]) = | | \hat{φ} - \hat{S} [f]) | |_{L_{2}},

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

\hat{φ}

$\hat{\varphi}$ - die beste Lösung. Dann werden die Verluste für unsere gewählte Strategie durch die Risikofunktion bestimmt :

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

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

Hier

P (f | φ)

$P(f|\varphi)$ bestimmt die Wahrscheinlichkeitsdichte unseres Ensembles, über die eine Verlustmittelung durchgeführt wird. Dieses Ensemble besteht aus einer hypothetischen Mehrfachwiederholung von Messungen.

f

$f$ für ein gegebenes

φ

$\varphi$ . Auf diese Weise,

P (f | φ)

$P(f|\varphi)$ - Dies ist die uns bekannte Wahrscheinlichkeitsdichte

f

$f$ im Experiment erhalten.

Nach dem Bayes'schen Ansatz wird vorgeschlagen, dies zu berücksichtigen

φ

$\varphi$ als Zufallsvariable mit a priori Wahrscheinlichkeitsdichte

P (φ)

$P(\varphi)$ Ausdruck der Zuverlässigkeit verschiedener Lösungen für unser Problem.

P (φ)

$P(\varphi)$ bestimmt auf der Grundlage von Informationen, die vor dem Experiment vorliegen. Dann basiert die Wahl der optimalen Strategie auf der Minimierung eines nachträglichen Risikos :

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

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

In diesem Fall ist die optimale Strategie bekannt:

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

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

Wo ist die hintere Dichte?

P (φ | f)

$P(\varphi|f)$ wird durch den Bayes-Satz bestimmt:

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

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

Dieser Ansatz ermöglicht es uns, die Varianz (Korrelationsfunktion) der resultierenden Lösung zu bestimmen:

D (x_{1}, x_{2}) = E [φ (x_{1}) - \hat{S} [f] (x_{1})] [φ (x_{2}) - \hat{S} [f] (x_{2})]

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

Wir haben also die optimale Lösung für unser Problem erhalten, indem wir eine a priori Dichte eingeführt haben

P (φ)

$P(\varphi)$ . Können wir etwas über diese Funktionswelt sagen?

φ (x)

$\varphi(x)$ was ist durch a priori Dichte gegeben?

Wenn die Antwort auf diese Frage Nein lautet, müssen wir alles Mögliche akzeptieren

φ (x)

$\varphi(x)$ ebenso wahrscheinlich und zu einer unregelmäßigen Lösung zurückkehren. Daher müssen wir diese Frage bejahen.

Genau darin besteht die Methode der statistischen Regularisierung - der Regularisierung einer Lösung durch Einführung zusätzlicher A-priori-Informationen zu

φ (x)

$\varphi(x)$ . Wenn der Forscher bereits a priori Informationen hat (a priori Dichte

P (\vec{φ})

$P(\vec{\varphi})$ ) kann er einfach das Integral berechnen und die Antwort bekommen.

Wenn es keine solchen Informationen gibt, beschreibt der nächste Absatz, welche minimalen Informationen ein Forscher möglicherweise hat und wie er sie verwendet, um eine regulierte Lösung zu erhalten.

A priori Informationen

Wie britische Wissenschaftler gezeigt haben, differenzieren sie im Rest der Welt gerne. Wenn der Mathematiker Fragen zur Rechtmäßigkeit dieser Operation stellt, glaubt der Physiker optimistisch, dass die Naturgesetze durch "gute" Funktionen beschrieben werden, dh glatt.

Mit anderen Worten, es macht es glatter

φ (x)

$\varphi(x)$ höhere a priori Wahrscheinlichkeitsdichte. Versuchen wir also, eine a priori-Wahrscheinlichkeit einzuführen, die auf der Glätte basiert. Zu diesem Zweck erinnern wir uns daran, dass die Einführung von A-priori-Informationen eine gewisse Gewalt gegen die Welt darstellt und die Naturgesetze dazu zwingt, auf eine für uns bequeme Weise auszusehen.

Diese Gewalt sollte minimiert werden, und durch die Einführung einer a priori Wahrscheinlichkeitsdichte ist es notwendig, dass Shannons Informationen über

φ (x)

$\varphi(x)$ Enthalten in

P (\vec{φ})

$P(\vec{\varphi})$ war minimal. Wenn wir das Obige formalisieren, leiten wir die Form der a priori-Dichte basierend auf der Glätte der Funktion ab. Dazu suchen wir nach einem bedingten Extrem an Informationen:

I [P (\vec{φ})] = \int \ln P (\vec{φ}) P (\vec{φ}) d \vec{φ} \to m i n

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

Unter folgenden Bedingungen:

Glättezustand $\varphi(x)$ . Lassen $\Omega$ Ist eine bestimmte Matrix, die die Glätte der Funktion kennzeichnet. Dann müssen wir einen bestimmten Wert der Glättefunktion erreichen:
$\int (\vec{φ}, Ω \vec{φ}) P (\vec{φ}) d \vec{φ} = ω$
$\int(\vec{\varphi},\Omega\vec{\varphi}) P(\vec{\varphi}) d\vec{\varphi} = \omega$
Ein aufmerksamer Leser sollte eine Frage zur Bestimmung des Werts eines Parameters stellen.
$\omega$ . Die Antwort wird im Text weiter gegeben.
Wahrscheinlichkeitsnormalisierung pro Einheit:
Unter diesen Bedingungen bietet die folgende Funktion ein Minimum an Funktionalität:
$P_{α} (\vec{φ}) = \frac{α^{R g (Ω) / 2} det Ω^{1 / 2}}{(2 π)^{N / 2}} \exp (- \frac{1}{2} (\vec{φ}, α Ω \vec{φ}))$
Parameter Verbunden mit Da wir jedoch keine Informationen über die spezifischen Werte der Glättungsfunktion haben, ist es sinnlos, genau herauszufinden, wie sie verbunden ist. Was ist dann damit zu tun?, du fragst. Hier werden Ihnen drei Wege offenbart:
1. Parameterwert anpassen $\alpha$ manuell und damit tatsächlich zu Tikhonovs Regularisierung übergehen
2. Mittelwertbildung (Integration) über alles Mögliche $\alpha$ unter der Annahme, dass alles möglich ist $\alpha$ ebenso wahrscheinlich
3. Wählen Sie die wahrscheinlichste $\alpha$ durch seine hintere Wahrscheinlichkeitsdichte $P(\alpha|\vec{f})$ . Dieser Ansatz ist korrekt, da er eine gute Annäherung an das Integral liefert, wenn die experimentellen Daten genügend Informationen über enthalten $\alpha$ .

Der erste Fall interessiert uns wenig. Im zweiten Fall müssen wir hier ein so hässliches Integral berechnen:

⟨ φ_{i} ⟩ = \frac{\int d φ φ_{i} P (f | φ) \int d α P (α) α^{\frac{R g (Ω)}{2}} \exp (- \frac{α}{2} (\vec{φ}, Ω \vec{φ}))}{\int d φ P (f | φ) \int d α P (α) α^{\frac{R g (Ω)}{2}} \exp (- \frac{α}{2} (\vec{φ}, Ω \vec{φ}))}

$\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}))}$

Für den dritten Fall können wir den Wert des Integrals analytisch für Gaußsche Rauschen im Experiment erhalten (dies wird im Abschnitt betrachtet).

Es sollte auch beachtet werden, dass wir das nirgendwo benutzt haben

Ω

$\Omega$ Ist ein Glättungsoperator. Tatsächlich können wir hier jeden anderen Operator (oder dessen lineare Kombination) verwenden. Nur die Glätte der Funktion ist die offensichtlichste Form von A-priori-Informationen, die wir verwenden können.

Probenahme

Wir haben über Funktionen gesprochen, aber jedes echte Gerät kann nicht nur ein Kontinuum messen, sondern auch eine zählbare Menge von Punkten. Wir nehmen immer Messungen in einer endlichen Menge von Punkten vor, daher sind wir gezwungen, Diskretisierungs- und Übergangsverfahren von der Integralgleichung zur Matrixgleichung durchzuführen. Bei der Methode der statistischen Regularisierung gehen wir wie folgt vor: Wir werden zerlegen

φ (x)

$\varphi(x)$ über ein System von Funktionen

{T_{n}}

$\{T_n\}$ ::

φ (x) = \sum_{n} φ_{n} T_{n} (x) .

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

Somit bilden die Koeffizienten dieser Expansion einen Vektor

\vec{φ}

$\vec{\varphi}$ Das ist ein Vektor im Funktionsraum.

Als funktionalen Raum können wir den Hilbert-Raum oder zum Beispiel den Raum der Polynome nehmen. Darüber hinaus ist die Wahl der Basis in diesen Räumen nur durch Ihre Vorstellungskraft begrenzt (wir haben versucht, mit den trigonometrischen Fourier-Reihen, Poligandra und kubischen Splines zu arbeiten).

Dann die Elemente der Matrix

K

$K$ berechnet als:

K_{m n} = (\hat{K} T_{n} (x)) (y_{m}),

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

y_{m}

$y_m$ - Punkte, an denen Messungen durchgeführt wurden. Matrixelemente

Ω

$\Omega$ wir berechnen nach der Formel:

Ω_{i j} = \int_{a}^{b} (\frac{d^{p} T_{i} (x)}{d x}) (\frac{d^{p} T_{j} (x)}{d x}) d x,

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

a

$a$ und

b

$b$ - die Grenzen des Intervalls, in dem die Funktion definiert ist

φ (x)

$\varphi(x)$ .

Verwenden Sie zur Neuberechnung von Fehlern die Dispersionsformel einer linearen Kombination von Zufallsvariablen:

D [φ (x)] = D [\sum_{n} φ_{n} T_{n} (x)] = \sum_{i, j} φ_{i} φ_{j} c o v (T_{i} (x), T_{j} (x)) .

$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)).$

Es ist zu beachten, dass in einigen Fällen die Darstellung einer Funktion unter Verwendung eines Vektors endlicher Dimension zu einem teilweisen Verlust oder einer Änderung von Informationen führt. Tatsächlich können wir Algebraisierung als eine Art Regularisierung betrachten, die jedoch schwach und unzureichend ist, um aus einer falschen Aufgabe eine richtige zu machen. Aber trotzdem haben wir uns jetzt von der Suche entfernt

φ (x)

$\varphi(x)$ zur Vektorsuche

\vec{φ}

$\vec{\varphi}$ und im nächsten Abschnitt finden wir es.

Gaußscher Rauschfall

Der Fall, in dem die Fehler im Experiment nach Gauß verteilt werden, ist insofern bemerkenswert
, als eine analytische Lösung für unser Problem erhalten werden kann. Da a priori Informationen und Fehler eine Gaußsche Form haben, hat ihr Produkt auch eine Gaußsche Form, und dann ist das hässliche Integral, das wir oben geschrieben haben, leicht zu nehmen. Die Lösung und ihr Fehler lauten wie folgt:

\vec{φ} = (K^{T} Σ^{- 1} K + α^{*} Ω)^{- 1} K^{T} Σ^{- 1^{T}} \vec{f}

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

Σ_{\vec{φ}} = (K^{T} Σ^{- 1} K + α^{*} Ω)^{- 1},

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

Σ

$\Sigma$ - Kovarianzmatrix der mehrdimensionalen Gaußschen Verteilung,

α^{*}

$\alpha^*$ - der wahrscheinlichste Wert des Parameters

α

$\alpha$ , die aus der Bedingung der maximalen a posteriori Wahrscheinlichkeitsdichte bestimmt wird:

P (α | \vec{f}) = C^{'} α^{\frac{R g (Ω)}{2}} \sqrt{| (K^{T} Σ^{- 1} K + α Ω)^{- 1} |} \exp (\frac{1}{2} {\vec{f}}^{T} Σ^{- 1} K^{T} (K^{T} Σ^{- 1} K + α Ω)^{- 1} K^{T} Σ^{- 1^{T}} \vec{f})

$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})$

Und wenn ich keine Gaußschen Fehler habe?

Der zweite Teil des Artikels wird diesem Thema gewidmet sein, aber lassen Sie uns zunächst das Wesentliche des Problems skizzieren.

⟨ φ_{i} ⟩ = \frac{\int d φ φ_{i} P (f | φ) \int d α P (α) α^{\frac{R g (Ω)}{2}} \exp (- \frac{α}{2} (\vec{φ}, Ω \vec{φ}))}{\int d φ P (f | φ) \int d α P (α) α^{\frac{R g (Ω)}{2}} \exp (- \frac{α}{2} (\vec{φ}, Ω \vec{φ}))}

Das Hauptproblem ist, dass dieses schreckliche Integral erstens mehrdimensional und zweitens in unendlichen Grenzen ist. Darüber hinaus ist es ein sehr mehrdimensionaler Vektor

\vec{φ}

$\vec{\varphi}$ kann leicht Dimension haben

m = 30 - 50

$m = 30-50$ und Gittermethoden zur Berechnung von Integralen weisen eine Komplexität des Typs auf

O (n^{m})

$O(n^m)$ daher in diesem Fall nicht anwendbar. Bei der Verwendung mehrdimensionaler Integrale funktioniert die Monte-Carlo-Integration gut.

Da unsere Grenzen unendlich sind, müssen wir außerdem die wichtige Stichprobenmethode verwenden, aber dann müssen wir eine Funktion für die Stichprobe auswählen. Um alles automatisierter zu machen, sollten Sie Markov Chain Monte Carlo (MCMC) verwenden , mit dem die Abtastfunktion unabhängig an das berechnete Integral angepasst werden kann. Wir werden im nächsten Artikel über die Anwendung von MCMC sprechen.

Praktischer Teil

Die erste Implementierung der statistischen Regularisierungsmethode wurde in den 70er Jahren auf Algol geschrieben und erfolgreich für Berechnungen in der atmosphärischen Physik verwendet. Trotz der Tatsache, dass wir immer noch die handschriftlichen Quellen des Algorithmus haben, haben wir uns entschlossen, ein bisschen Modernismus hinzuzufügen und eine Implementierung in Python und dann in Julia vorzunehmen.

Python

Installation

Installieren über pip:

pip install statreg

oder laden Sie die Quelle Code .

Beispiele

Betrachten Sie als Beispiel, wie Sie mit einem Modul staregDaten für eine Matrix und eine Integralgleichung wiederherstellen.

Wir importieren die notwendigen wissenschaftlichen Pakete.

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

Wir bestimmen das wahre Signal, das wir wiederherstellen werden.

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));

Definieren Sie den Kernel und die Funktionsweise der Faltung von Funktionen (Hinweis: np.convolutionspeziell für Arrays):

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])

Wir erzeugen die gemessenen Daten und rauschen sie unter Verwendung der Normalverteilung:

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);

Wir lösen die Integralgleichung

Wir importieren den Solver und die Hilfsklasse zur Diskretisierung:

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

Als funktionale Basis für die Diskretisierung verwenden wir kubische Splines und als zusätzliche Bedingung geben wir an, dass die Funktion an den Kanten Nullwerte annimmt.

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

Wir lösen die Gleichung:

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

Wir erstellen einen Lösungsplan:

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);

Wir lösen die Matrixgleichung

Wir importieren den Solver und die Hilfsklasse zur Diskretisierung:

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

Um Matrizen zu erhalten, verwenden wir unsere funktionale Basis, aber es ist klar, dass Matrizen auf jede Weise erhalten werden können.

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

Wir lösen die Matrixgleichung:

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

Erstellen Sie das Diagramm:

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

Wie bereits erwähnt, erfordert die Weiterentwicklung der Technik eine erweiterte Monte-Carlo-Integration. Wir könnten ein Modul in Python verwenden (zum Beispiel haben wir mit PyMC3 gearbeitet), aber wir beteiligen uns unter anderem an einem gemeinsamen Projekt mit dem Max-Planck-Institut in München.

Dieses Projekt heißt Bayesian Analysis Toolkit . Ziel ist es, ein Framework mit Tools für Bayes'sche Analysemethoden zu schaffen, einschließlich Tools für MCMC. Jetzt arbeitet das Team an der zweiten Version des Frameworks, die in Julia geschrieben ist (die erste ist in schlechtem C ++ geschrieben). Eine der Aufgaben unserer Gruppe ist es, die Fähigkeiten dieses Frameworks am Beispiel der statistischen Regularisierung zu demonstrieren. Deshalb haben wir eine Implementierung in Julia geschrieben .

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)

Wenn wir diesmal einen anderen Kern verwenden, werden wir keinen Integrationsschritt machen, sondern eine Faltung mit einem Gaußschen, was tatsächlich zu einer gewissen „Unschärfe“ unserer Daten führt:

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)

Ebenso nehmen wir die Basis von Splines mit festen Enden:

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)

Beispiel aus der Praxis

Als Beispiel für die Analyse realer Daten werden wir das Elektronenstreuspektrum eines Wasserstoff-Deuterium-Gemisches wiederherstellen. Im Experiment wurde das Integralspektrum gemessen (dh die Anzahl der Elektronen liegt über einer bestimmten Energie), und wir müssen das Differenzspektrum wiederherstellen. Für diese Daten wurde das Spektrum zunächst mithilfe der Anpassung rekonstruiert, sodass wir eine Grundlage für die Überprüfung der Richtigkeit unseres Algorithmus haben.

So sieht das anfänglich integrierte Spektrum aus:

Und so - das Ergebnis der Restauration: Die

Analyse mit Anpassung hat drei Hauptnachteile:

Die statistische Regularisierung vermeidet all diese Probleme und liefert ein modellunabhängiges Ergebnis mit Messfehlern. Die durch Regularisierung erhaltene Lösung stimmt gut mit der Anpassungskurve überein. Beachten Sie die beiden kleinen Peaks bei 25 und 30 eV. Es ist bekannt, dass während der Doppelstreuung ein Peak bei 25 eV gebildet wird, der durch eine Anpassung wiederhergestellt wurde, da er in der Anpassungsfunktion eindeutig spezifiziert wurde. Ein Peak von 30 eV kann eine statistische Anomalie sein (die Fehler sind zu diesem Zeitpunkt ziemlich groß) oder möglicherweise auf das Vorhandensein zusätzlicher dissoziativer Streuung hinweisen.

Schlussfolgerungen und Ankündigung des nächsten Teils

Wir haben Ihnen eine nützliche Technik vorgestellt, die an viele Aufgaben der Datenanalyse (einschließlich maschinelles Lernen) angepasst werden kann und eine ehrliche "Übereinstimmung" der Antwort liefert - die rationalste Lösung für die Gleichung angesichts der durch Messfehler verursachten Unsicherheit. Als netten Bonus erhalten wir Werte für den Entscheidungsfehler. Wer an der Entwicklung teilnehmen oder die Methode der statistischen Regularisierung anwenden möchte, kann in Form von Code in Python, Julia oder auf etwas anderem beitragen.

Im nächsten Teil werden wir sprechen über:

MCMC verwenden
Cholesky-Zersetzung
Als praktisches Beispiel betrachten wir die Verwendung der Regularisierung zur Verarbeitung eines Signals aus einem Modell eines Orbitaldetektors von Protonen und Elektronen

Verweise

Gepostet von Mikhail Zeleny , Forscher am Labor für Kernphysik-Versuchsmethoden bei JetBrains Research .

Statistische Regularisierung falscher inverser Probleme zu ihnen. Turchin (Teil 1)