Régularisation statistique de problèmes inverses incorrects. Turchin (partie 1)

Bonjour, Habr! Aujourd'hui, nous voulons vous dire ce que fait le laboratoire de méthodes expérimentales de physique nucléaire , qui fait partie de JetBrains Research .

Où sont JetBrains et où est la physique nucléaire, demandez-vous. Nous nous sommes entendus sur la base de l'amour pour Kotlin, bien que dans ce post, nous ne parlerons pas de lui. Notre groupe se concentre sur le développement de logiciels d'analyse, de modélisation et d'écriture de données pour les scientifiques, et se concentre donc sur la coopération et le partage des connaissances avec les entreprises informatiques.

Dans cet article, nous voulons parler de la méthode de régularisation statistique que nous popularisons , proposée par V.F. Turchin dans les années 70 du XXe siècle, et de sa mise en œuvre sous forme de code en Python et Julia.

La présentation sera assez détaillée, donc ceux qui sont clairs sur les problèmes inverses peuvent aller directement aux exemples et lire la théorie de cet article .

• Occurrence du problème: pourquoi devrait-on régulariser du tout?
• Description théorique de la régularisation statistique
• Exemples pratiques sur:
  
  • Python
  • Julia
• Conclusion

Occurrence du problème: pourquoi devrait-on régulariser du tout?

S'il suffit d'ignorer, toute mesure de l'expérience peut être décrite comme suit: il existe un certain appareil qui capture le spectre ou le signal d'un processus et affiche certains nombres en fonction des résultats de la mesure. Notre tâche en tant que chercheurs, en examinant ces chiffres et en connaissant la structure de l'appareil, est de comprendre quel était le spectre ou le signal mesuré. C'est-à-dire face à ce qu'on appelle le problème inverse . Si vous l'imaginez mathématiquement, nous obtenons cette équation (qui, soit dit en passant , est appelée l'équation de Fredholm du premier type ):

$$

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

En fait, cette équation décrit ce qui suit: notre appareil de mesure est représenté ici par sa fonction matérielle $$$K(x,y)$qui agit sur le spectre étudié ou tout autre signal d'entrée $$$\varphi$obligeant le chercheur à observer le signal de sortie $$$f(y)$. L'objectif du chercheur est de restaurer le signal$$$\varphi$ par célèbre $$$f(y)$et $$$K(x,y)$. Vous pouvez également formuler cette expression sous forme matricielle, en remplaçant les fonctions par des vecteurs et des matrices:

$$

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

Il semblerait que la reconstruction du signal ne soit pas une tâche difficile, car l'équation de Fredholm et le système d'équations linéaires (même surdéterminés) ont une solution exacte. Alors essayons. Que le signal mesuré soit décrit comme la somme de deux gausses:

$$

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

En tant qu'instrument, nous prenons l'intégrateur le plus simple - une matrice qui traduit notre signal en une somme cumulative en utilisant la fonction Heaviside:

$$

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

Le type du signal mesuré et de notre appareil, ainsi que le résultat de la mesure sont indiqués sur le graphique.

Il est important que toute mesure réelle ait une erreur, donc nous gâcherons un peu notre résultat en ajoutant un bruit normal, donnant une erreur de mesure de cinq pour cent.

Nous allons restaurer le signal par la méthode des moindres carrés:

$$

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

Et en conséquence, nous obtenons:

En fait, sur ce point, nous pourrions terminer l'article, une fois de plus, nous étant convaincus de l'impuissance des méthodes mathématiques idéalistes face à la réalité physique dure et impitoyable, et aller fumer des fers à souder.

Mais d'abord, essayons de comprendre ce qui nous a causé cet échec? Évidemment, il s'agit d'erreurs de mesure, mais qu'est-ce qu'elles affectent? Le fait est que Jacques Hadamard (le même qui a ajouté un tiret à la formule Cauchy - Hadamard) a divisé toutes les tâches en tâches correctement posées et incorrectes.

Se souvenir des classiques: «Cela n'a aucun sens de chercher une solution, le cas échéant. Il s'agit de savoir comment gérer une tâche qui n'a pas de solution. C'est une question profondément fondamentale ... »- nous ne parlerons pas des tâches correctes et ne reprendrons pas immédiatement les tâches incorrectes. Heureusement, nous l'avons déjà rencontré: l'équation de Fredholm écrite ci-dessus est un problème inverse incorrect - même avec des fluctuations infiniment petites dans les données d'entrée (et même nos erreurs de mesure sont loin d'être infinitésimales), la solution de l'équation obtenue de manière analytique exacte peut différer arbitrairement de la vraie .

Vous pouvez lire la preuve de cette affirmation dans le premier chapitre du travail classique de l'académicien A.N. Tikhonova "Méthodes pour résoudre les problèmes mal posés." Ce livre contient des conseils sur ce qu'il faut faire avec des tâches incorrectes, mais la technique décrite présente un certain nombre d'inconvénients qui sont résolus dans la méthode Turchin. Mais d'abord, nous décrivons les principes généraux du travail avec des tâches incorrectes: que faire si vous rencontrez une telle tâche?

Étant donné que la tâche elle-même ne peut rien nous offrir, nous devrons commettre un petit délit: compléter la tâche avec des données afin qu'elle devienne correcte, en d'autres termes, saisir des * informations a priori supplémentaires sur la tâche * (ce processus est appelé régularisation de la tâche). Contrairement à la méthode classique de Tikhonov, basée sur l'introduction de fonctionnelles de régularisation paramétrisées, la méthode Turchin fait appel en revanche à des méthodes bayésiennes.

Description théorique de la régularisation statistique

Stratégie

Nous formulons notre problème en termes de statistiques mathématiques: selon une implémentation bien connue $$$f$(que nous mesurons dans l'expérience), nous devons évaluer la valeur du paramètre $$$\varphi$. Fonctionnel$$$\hat{S}$ calculateur $$$\varphi$basé $$$f$nous appellerons stratégie . Afin de déterminer les stratégies les plus optimales, nous introduisons une fonction de perte quadratique . La fonction de perte réelle peut être quelconque, pourquoi choisissons-nous celle quadratique? Parce que toute fonction de perte proche de son minimum peut être approximée par une fonction quadratique:

$$

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

Où $$$\hat{\varphi}$- la meilleure solution. Ensuite, les pertes pour notre stratégie choisie sont déterminées par la fonction de risque :

$$

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

Ici $$$P(f|\varphi)$détermine la densité de probabilité de notre ensemble, sur laquelle la moyenne des pertes est effectuée. Cet ensemble est formé par une répétition multiple hypothétique de mesures.$$$f$ pour une donnée $$$\varphi$. De cette façon,$$$P(f|\varphi)$ - c'est la densité de probabilité que nous connaissons $$$f$obtenus dans l'expérience.

Selon l'approche bayésienne, il est proposé de considérer$$$\varphi$comme une variable aléatoire avec une densité de probabilité a priori $$$P(\varphi)$exprimant la fiabilité de diverses solutions à notre problème.$$$P(\varphi)$déterminé sur la base des informations existantes avant l'expérience. Ensuite, le choix de la stratégie optimale est basé sur la minimisation du risque a posteriori :

$$

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

Dans ce cas, la stratégie optimale est bien connue:

$$

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

où est la densité postérieure $$$P(\varphi|f)$est déterminé par le théorème de Bayes:

$$

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

Cette approche nous permettra de déterminer la variance (fonction de corrélation) de la solution résultante:

$$

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

Nous avons donc obtenu la solution optimale à notre problème en introduisant une densité a priori $$$P(\varphi)$. Pouvons-nous dire quoi que ce soit sur ce monde de fonctions$$$\varphi(x)$qui est donné par la densité a priori?

Si la réponse à cette question est non, nous devrons accepter tous les$$$\varphi(x)$tout aussi probable et retour à une solution irrégulière. Par conséquent, nous devons répondre à cette question par l'affirmative.

C’est précisément ce que consiste la méthode de régularisation statistique - régularisation d’une solution en introduisant des informations a priori$$$\varphi(x)$. Si le chercheur dispose déjà d'informations a priori (densité a priori$$$P(\vec{\varphi})$), il peut simplement calculer l'intégrale et obtenir la réponse.

S'il n'y a pas de telles informations, le paragraphe suivant décrit les informations minimales qu'un chercheur peut avoir et comment les utiliser pour obtenir une solution régularisée.

Information a priori

Comme l'ont montré les scientifiques britanniques, dans le reste du monde, ils aiment se différencier. De plus, si le mathématicien pose des questions sur la légalité de cette opération, le physicien estime avec optimisme que les lois de la nature sont décrites par de «bonnes» fonctions, c'est-à-dire lisses.

En d'autres termes, cela rend plus fluide$$$\varphi(x)$densité de probabilité a priori plus élevée. Essayons donc d'introduire une probabilité a priori basée sur la fluidité. Pour ce faire, nous rappelons que l'introduction de l'information a priori est une violence contre le monde, forçant les lois de la nature à nous paraître d'une manière qui nous convient.

Cette violence doit être minimisée et, en introduisant une densité de probabilité a priori, il est nécessaire que les informations de Shannon sur$$$\varphi(x)$contenu dans $$$P(\vec{\varphi})$était minime. En formalisant ce qui précède, nous dérivons la forme d'une densité a priori basée sur la finesse de la fonction. Pour ce faire, nous chercherons un extremum conditionnel d'information:

$$

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

Dans les conditions suivantes:

1. Condition de douceur $$$\varphi(x)$. Laisser être$$$\Omega$Est une certaine matrice caractérisant la fluidité de la fonction. Ensuite, nous exigeons qu'une certaine valeur de la fonction de lissage soit atteinte:
   

   $$

   $$\int(\vec{\varphi},\Omega\vec{\varphi}) P(\vec{\varphi}) d\vec{\varphi} = \omega$$
   Un lecteur attentif devrait poser une question sur la détermination de la valeur d'un paramètre.
   $$$\omega$. La réponse sera donnée plus loin dans le texte.
   
2. Normalisation des probabilités par unité: $$$\int P(\vec{\varphi}) d\vec{\varphi} = 1$
   Dans ces conditions, la fonction suivante fournira un minimum de fonctionnalités:
   

   $$

   $$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}))$$
   Paramètre $$$\alpha$Connecté avec $$$\omega$, mais comme nous ne disposons pas d'informations sur les valeurs spécifiques de la fonction de lissage, il est inutile de savoir exactement comment elle est connectée. Que faire ensuite avec$$$\alpha$, tu demandes. Voici trois façons de vous révéler:
   
   
   1. Correspond à la valeur du paramètre $$$\alpha$manuellement et ainsi procéder à la régularisation de Tikhonov
   2. Moyennage (intégration) sur tous les possibles $$$\alpha$en supposant tout possible $$$\alpha$tout aussi probable
      
   3. Choisissez le plus probable $$$\alpha$par sa densité de probabilité postérieure $$$P(\alpha|\vec{f})$. Cette approche est correcte car elle donne une bonne approximation de l'intégrale si les données expérimentales contiennent suffisamment d'informations sur$$$\alpha$.

Le premier cas nous intéresse peu. Dans le second cas, nous devons calculer une telle intégrale laide ici:

$$

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

Pour le troisième cas, nous pouvons obtenir la valeur de l'intégrale analytiquement pour les bruits gaussiens dans l'expérience (cela sera considéré dans la section).

Il convient également de noter que nous n'avons utilisé nulle part, que$$$\Omega$Est un opérateur de douceur. En fait, nous pouvons utiliser n'importe quel autre opérateur (ou leur combinaison linéaire) ici, juste la finesse de la fonction est la forme la plus évidente d'informations a priori que nous pouvons utiliser.

Échantillonnage

Nous avons parlé de fonctions, mais tout appareil réel ne peut pas mesurer non seulement un continuum, mais même un ensemble dénombrable de points. Nous prenons toujours des mesures dans un ensemble fini de points, nous sommes donc obligés d'effectuer des procédures de discrétisation et de transition de l'équation intégrale à celle de la matrice. Dans la méthode de régularisation statistique, nous procédons comme suit: nous décomposerons$$$\varphi(x)$ sur un certain système de fonctions $$$\{T_n\}$:

$$

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

Ainsi, les coefficients de cette expansion forment un vecteur $$$\vec{\varphi}$qui est un vecteur dans l'espace des fonctions.

En tant qu'espace fonctionnel, nous pouvons prendre l'espace de Hilbert, ou, par exemple, l'espace des polynômes. De plus, le choix de la base dans ces espaces n'est limité que par votre imagination (nous avons essayé de travailler avec la série trigonométrique de Fourier, la poligandre et les splines cubiques).

Ensuite, les éléments de la matrice$$$K$ calculé comme:

$$

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

Où $$$y_m$- points auxquels les mesures ont été effectuées. Éléments de matrice$$$\Omega$ nous calculerons par la formule:

$$

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

Où $$$a$et $$$b$- les limites de l'intervalle sur lequel la fonction est définie $$$\varphi(x)$.

Pour recalculer les erreurs, utilisez la formule de dispersion d'une combinaison linéaire de variables aléatoires:

$$

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

Il faut garder à l'esprit que dans certains cas, la représentation d'une fonction à l'aide d'un vecteur de dimension finie entraîne une perte ou un changement partiel d'information. En fait, nous pouvons considérer l'algèbre comme une sorte de régularisation, même faible et insuffisante pour transformer une tâche incorrecte en une tâche correcte. Mais, de toute façon, nous sommes maintenant passés de la recherche$$$\varphi(x)$ recherche vectorielle $$$\vec{\varphi}$et dans la section suivante, nous le trouvons.

Cas de bruit gaussien

Le cas où les erreurs de l'expérience sont réparties selon Gauss est remarquable
en ce qu'une solution analytique à notre problème peut être obtenue. Puisque les informations et les erreurs a priori ont une forme gaussienne, leur produit a également une forme gaussienne, puis l'intégrale laide que nous avons écrite ci-dessus est facilement prise. La solution et son erreur seront les suivantes:

$$

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

Où $$$\Sigma$- matrice de covariance de la distribution gaussienne multidimensionnelle, $$$\alpha^*$- la valeur la plus probable du paramètre $$$\alpha$, qui est déterminée à partir de la condition de densité de probabilité maximale 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})$$

Et si je n'ai pas d'erreurs gaussiennes?

La deuxième partie de l'article sera consacrée à cela, mais pour l'instant résumons l'essence du problème.

$$

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

Le problème principal est que cette terrible intégrale, d'une part multidimensionnelle, et d'autre part, dans des limites infinies. De plus, il est hautement multidimensionnel, vecteur$$$\vec{\varphi}$ peut facilement avoir une dimension $$$m = 30-50$et les méthodes de grille pour calculer les intégrales ont une complexité de type $$$O(n^m)$, par conséquent, inapplicable en l'espèce. Lors de la prise d'intégrales multidimensionnelles, l'intégration Monte Carlo fonctionne bien.

De plus, comme nos limites sont infinies, nous devons utiliser la méthode d'échantillonnage importante, mais nous devons ensuite sélectionner une fonction d'échantillonnage. Pour rendre tout plus automatisé, vous devez utiliser la chaîne de Markov Monte Carlo (MCMC) , qui peut adapter indépendamment la fonction d'échantillonnage à l'intégrale calculée. Nous parlerons de l'application de MCMC dans le prochain article.

Partie pratique

La première mise en œuvre de la méthode de régularisation statistique a été écrite dans les années 70 sur Algol et a été utilisée avec succès pour des calculs en physique atmosphérique. Malgré le fait que nous ayons encore les sources manuscrites de l'algorithme, nous avons décidé d'ajouter un peu de modernisme et de faire une implémentation en Python, puis sur Julia.

Python

Installation

Installer via pip:

pip install statreg

ou téléchargez le code source .

Exemples

Par exemple, considérez comment utiliser un module staregpour récupérer des données pour une matrice et une équation intégrale.

Nous importons les packages scientifiques nécessaires.

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

Nous déterminons le vrai signal, que nous rétablirons.

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

Définir le noyau et le fonctionnement de la convolution des fonctions (Remarque: np.convolutionspécifiquement pour les tableaux):

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

Nous générons les données mesurées et les faisons du bruit en utilisant la distribution normale:

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

Nous résolvons l'équation intégrale

Nous importons le solveur et la classe auxiliaire pour la discrétisation:

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

Comme base fonctionnelle pour la discrétisation, nous utilisons des splines cubiques, et comme condition supplémentaire, nous indiquons que la fonction prend des valeurs nulles sur les bords.

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

Nous résolvons l'équation:

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

Nous construisons un calendrier de solutions:

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

Nous résolvons l'équation de la matrice

Nous importons le solveur et la classe auxiliaire pour la discrétisation:

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

Pour obtenir des matrices, nous utilisons notre base fonctionnelle, mais il est clair que les matrices peuvent être obtenues de n'importe quelle manière.

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

Nous résolvons l'équation matricielle:

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

Construisez le graphique:

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

Comme nous l'avons mentionné, la poursuite du développement de la technique nécessite une intégration Monte Carlo avancée. Nous pourrions utiliser un module en Python (par exemple, nous avons travaillé avec PyMC3), cependant, nous participons, entre autres, à un projet conjoint avec le Max Planck Institute de Munich.

Ce projet s'appelle le Bayesian Analysis Toolkit . Son objectif est de créer un cadre avec des outils pour les méthodes d'analyse bayésienne, comprenant principalement des outils pour MCMC. Maintenant, l'équipe travaille sur la deuxième version du framework, qui est écrite en Julia (la première est écrite en mauvais C ++). L'une des tâches de notre groupe est de démontrer les capacités de ce cadre en utilisant l'exemple de la régularisation statistique, nous avons donc écrit une implémentation dans 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)

Cette fois, nous utilisons un noyau différent, nous ne prendrons pas une étape d'intégration, mais une convolution avec un gaussien, ce qui conduit en fait un certain «flou» à nos données:

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)

De même, nous prenons la base de splines à extrémités fixes:

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)

Exemple du monde réel

À titre d'exemple d'analyse de données réelles, nous restaurerons le spectre de diffusion d'électrons d'un mélange hydrogène-deutérium. Dans l'expérience, le spectre intégral a été mesuré (c'est-à-dire que le nombre d'électrons est supérieur à une certaine énergie), et nous devons restaurer le spectre différentiel. Pour ces données, le spectre a été initialement reconstruit en utilisant l'ajustement, de sorte que nous avons une base pour vérifier l'exactitude de notre algorithme.

Voici à quoi ressemble le spectre intégré initial:

Et donc - le résultat de la restauration:

L'analyse avec ajustement présente trois principaux inconvénients:

• , .
• , .
• .

La régularisation statistique évite tous ces problèmes et fournit un résultat indépendant du modèle avec des erreurs de mesure. La solution obtenue par régularisation est en bon accord avec la courbe d'ajustement. Notez les deux petits pics à 25 et 30 eV. Il est connu qu'un pic à 25 eV est formé lors de la double diffusion, et il a été restauré par un ajustement, car il était clairement spécifié dans la fonction d'ajustement. Un pic de 30 eV peut être une anomalie statistique (les erreurs sont assez importantes à ce stade), ou, éventuellement, indique la présence d'une diffusion dissociative supplémentaire.

Conclusions et annonce de la prochaine partie

Nous vous avons parlé d'une technique utile qui peut être adaptée à de nombreuses tâches d'analyse de données (y compris l'apprentissage automatique), et obtenir un "ajustement" honnête de la réponse - la solution la plus rationnelle à l'équation face à l'incertitude causée par les erreurs de mesure. En prime, nous obtenons des valeurs pour l'erreur de décision. Ceux qui veulent participer au développement ou appliquer la méthode de régularisation statistique peuvent contribuer sous forme de code en Python, Julia ou sur autre chose.

Dans la partie suivante, nous parlerons de:

• Utilisation de MCMC
• Décomposition Cholesky
• À titre d'exemple pratique, nous considérons l'utilisation de la régularisation pour le traitement d'un signal provenant d'un modèle de détecteur orbital de protons et d'électrons
  

Références

• Oeuvre originale Turchin
• La thèse sur l'application de la régularisation en physique atmosphérique
• Notre publication sur l'application de la régularisation statistique
• Implémentation de Python
• Implémentation sur Julia
• Implémentation sur ALGOL-60 , pour les amateurs d'archéologie
• Boîte à outils d'analyse bayésienne

Publié par Mikhail Zeleny , chercheur au laboratoire de méthodes d'expérimentation en physique nucléaire de JetBrains Research .