Regularização estatística de problemas inversos incorretos para eles. Turchin (parte 1)

Olá Habr! Hoje, queremos contar a você o que o laboratório de métodos de experimentos de física nuclear , parte da JetBrains Research , faz .

Onde estão os JetBrains e onde está a física nuclear, você pergunta. Concordamos com base no amor por Kotlin, embora neste post não falaremos sobre ele. Nosso grupo se concentra no desenvolvimento de software de análise de dados, modelagem e gravação para cientistas e, portanto, se concentra na cooperação e no compartilhamento de conhecimento com empresas de TI.

Neste artigo, queremos falar sobre o método de regularização estatística que popularizamos , proposto por V.F. Turchin nos anos 70 do século XX, e sua implementação na forma de código em Python e Julia.

A apresentação será bastante detalhada, para que aqueles que são claros sobre os problemas inversos possam ir direto aos exemplos e ler a teoria neste artigo .

• Ocorrência do problema: por que alguém deveria se regularizar?
• Descrição Teórica da Regularização Estatística
• Exemplos práticos sobre:
  
  • Pitão
  • Julia
• Conclusão

Ocorrência do problema: por que alguém deveria se regularizar?

Se for suficiente ignorar, qualquer medição no experimento pode ser descrita da seguinte forma: existe um determinado dispositivo que captura o espectro ou sinal de um processo e mostra alguns números de acordo com os resultados da medição. Nossa tarefa como pesquisadores, observando esses números e conhecendo a estrutura do dispositivo, é entender qual era o espectro ou sinal medido. Ou seja, diante do que é chamado de problema inverso . Se você imaginar isso matematicamente, obteremos esta equação (que, aliás, é chamada de equação de Fredholm do primeiro tipo ):

$$

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

De fato, esta equação descreve o seguinte: nosso dispositivo de medição é representado aqui por sua função de hardware $$$K(x,y)$que atua no espectro estudado ou em outro sinal de entrada $$$\varphi$fazendo com que o pesquisador observe o sinal de saída $$$f(y)$. O objetivo do pesquisador é restaurar o sinal$$$\varphi$ pelo famoso $$$f(y)$e $$$K(x,y)$. Você também pode formular esta expressão na forma de matriz, substituindo funções por vetores e matrizes:

$$

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

Parece que a reconstrução do sinal não é uma tarefa difícil, pois tanto a equação de Fredholm quanto o sistema de equações lineares (mesmo superdeterminadas) têm uma solução exata. Então, vamos tentar. Deixe o sinal medido ser descrito como a soma de dois gausses:

$$

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

Como instrumento, tomamos o integrador mais simples - uma matriz que traduz nosso sinal em uma soma cumulativa usando a função Heaviside:

$$

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

O tipo do sinal medido e nosso dispositivo, bem como o resultado da medição são mostrados no gráfico.

É importante que qualquer medição real tenha um erro, portanto, estragaremos um pouco o resultado adicionando ruído normal, resultando em um erro de medição de cinco por cento.

Restauraremos o sinal pelo método dos mínimos quadrados:

$$

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

E, como resultado, obtemos: na

verdade, poderíamos terminar o artigo, mais uma vez convencidos do desamparo dos métodos idealistas da matemática diante da realidade física dura e cruel, e começamos a fumar ferros de solda.

Mas primeiro, vamos descobrir o que causou esse fracasso conosco. Obviamente, o ponto são os erros de medição, mas o que eles afetam? O fato é que Jacques Hadamard (o mesmo que adicionou uma pitada à fórmula de Cauchy - Hadamard) dividiu todas as tarefas em tarefas corretamente colocadas e incorretas.

Lembrando os clássicos: “Não faz sentido procurar uma solução, se houver. É sobre como lidar com uma tarefa que não tem solução. Esta é uma questão profundamente fundamental ... ”- não falaremos sobre as tarefas corretas e imediatamente assumiremos as incorretas. Felizmente, já encontramos o seguinte: a equação de Fredholm escrita acima é um problema inverso incorreto - mesmo com flutuações infinitamente pequenas nos dados de entrada (e até nossos erros de medição estão longe de serem infinitesimais), a solução da equação obtida de uma maneira analítica exata pode diferir arbitrariamente da verdadeira .

Você pode ler a prova dessa afirmação no primeiro capítulo da obra clássica do acadêmico A.N. Tikhonova "Métodos para resolver problemas incorretos". Este livro tem dicas sobre o que fazer com tarefas incorretas, mas a técnica descrita lá tem várias desvantagens que são corrigidas no método Turchin. Mas primeiro, descrevemos os princípios gerais do trabalho com tarefas incorretas: o que fazer se você se deparar com uma tarefa dessas?

Como a tarefa em si não pode oferecer nada, teremos que cometer um pequeno crime: completar a tarefa com dados para que ela se torne correta, ou seja, digite * informações adicionais a priori adicionais sobre a tarefa * (esse processo é chamado de regularização da tarefa). Em contraste com o método clássico de Tikhonov, baseado na introdução de funcionais de regularização parametrizados, o método de Turchin chama, por outro lado, usando métodos bayesianos.

Descrição Teórica da Regularização Estatística

Estratégia

Formulamos nosso problema em termos de estatística matemática: de acordo com uma implementação bem conhecida $$$f$(que medimos no experimento), precisamos avaliar o valor do parâmetro $$$\varphi$. Funcional$$$\hat{S}$ calculando $$$\varphi$Sediada $$$f$chamaremos de estratégia . Para determinar quais estratégias são mais ideais, introduzimos uma função de perda quadrática . A função de perda real pode ser qualquer uma, por que escolhemos a função quadrática? Como qualquer função de perda próxima ao mínimo pode ser aproximada por uma função quadrática:

$$

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

Onde $$$\hat{\varphi}$- a melhor solução. Em seguida, as perdas para a estratégia escolhida são determinadas pela função de risco :

$$

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

Aqui $$$P(f|\varphi)$determina a densidade de probabilidade de nosso conjunto, sobre a qual a média da perda é realizada. Esse conjunto é formado por uma repetição múltipla hipotética de medições.$$$f$ para um dado $$$\varphi$. Nesse caminho,$$$P(f|\varphi)$ - esta é a densidade de probabilidade conhecida por nós $$$f$obtido no experimento.

De acordo com a abordagem bayesiana, propõe-se considerar$$$\varphi$como uma variável aleatória com densidade de probabilidade a priori $$$P(\varphi)$expressando a confiabilidade de várias soluções para o nosso problema.$$$P(\varphi)$determinado com base nas informações existentes antes da experiência. Então a escolha da estratégia ideal é baseada na minimização do risco a posteriori :

$$

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

Nesse caso, a estratégia ideal é bem conhecida:

$$

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

onde é a densidade posterior $$$P(\varphi|f)$é determinado pelo teorema de Bayes:

$$

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

Essa abordagem nos permitirá determinar a variação (função de correlação) da solução resultante:

$$

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

Portanto, obtivemos a solução ideal para o nosso problema, introduzindo uma densidade a priori $$$P(\varphi)$. Podemos dizer algo sobre esse mundo de funções$$$\varphi(x)$qual é dado por uma densidade a priori?

Se a resposta a esta pergunta for não, teremos de aceitar todos os possíveis$$$\varphi(x)$igualmente provável e retornar a uma solução irregular. Portanto, devemos responder afirmativamente a essa pergunta.

É precisamente nisso que consiste o método de regularização estatística - regularização de uma solução, introduzindo informações a priori adicionais sobre$$$\varphi(x)$. Se o pesquisador já tiver alguma informação a priori (densidade a priori$$$P(\vec{\varphi})$), ele pode apenas calcular a integral e obter a resposta.

Se não houver essas informações, o próximo parágrafo descreve quais informações mínimas um pesquisador pode ter e como usá-las para obter uma solução regularizada.

Informação a priori

Como os cientistas britânicos demonstraram, no resto do mundo eles gostam de se diferenciar. Além disso, se o matemático fizer perguntas sobre a legalidade dessa operação, o físico acredita otimista que as leis da natureza são descritas por funções "boas", ou seja, suaves.

Em outras palavras, torna mais suave$$$\varphi(x)$maior densidade de probabilidade a priori. Então, vamos tentar introduzir uma probabilidade a priori com base na suavidade. Para fazer isso, lembramos que a introdução de informações a priori representa alguma violência contra o mundo, forçando as leis da natureza a parecerem de uma maneira conveniente para nós.

Essa violência deve ser minimizada e, ao introduzir uma densidade de probabilidade a priori, é necessário que as informações de Shannon sobre$$$\varphi(x)$contido em $$$P(\vec{\varphi})$foi mínimo. Formalizando o exposto, derivamos a forma de densidade a priori com base na suavidade da função. Para fazer isso, procuraremos um extremo condicional da informação:

$$

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

Sob as seguintes condições:

1. Condição de suavidade $$$\varphi(x)$. Deixe ser$$$\Omega$É uma certa matriz que caracteriza a suavidade da função. Então exigimos que um certo valor da suavidade funcional seja alcançado:
   

   $$

   $$\int(\vec{\varphi},\Omega\vec{\varphi}) P(\vec{\varphi}) d\vec{\varphi} = \omega$$
   Um leitor atento deve fazer uma pergunta sobre como determinar o valor de um parâmetro.
   $$$\omega$. A resposta será dada mais adiante no texto.
   
2. Normalização de probabilidade por unidade: $$$\int P(\vec{\varphi}) d\vec{\varphi} = 1$
   Sob essas condições, a seguinte função fornecerá um mínimo de funcionalidade:
   

   $$

   $$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}))$$
   Parâmetro $$$\alpha$Conectado com $$$\omega$, mas como não temos informações sobre os valores específicos da suavidade funcional, descobrir exatamente como ele está conectado é inútil. O que fazer então com$$$\alpha$, você pergunta. Aqui, três maneiras são reveladas a você:
   
   
   1. Corresponder ao valor do parâmetro $$$\alpha$manualmente e, assim, efetivamente proceder à regularização de Tikhonov
   2. Média (integração) de todos os possíveis $$$\alpha$assumindo tudo possível $$$\alpha$igualmente provável
      
   3. Escolha o mais provável $$$\alpha$por sua densidade de probabilidade posterior $$$P(\alpha|\vec{f})$. Essa abordagem está correta porque fornece uma boa aproximação da integral se os dados experimentais contiverem informações suficientes sobre$$$\alpha$.

O primeiro caso é de pouco interesse para nós. No segundo caso, devemos calcular uma integral tão feia aqui:

$$

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

Para o terceiro caso, podemos obter o valor da integral analiticamente para ruídos gaussianos no experimento (isso será considerado na seção).

Note-se também que não usamos em nenhum lugar, que$$$\Omega$É um operador de suavidade. De fato, podemos usar qualquer outro operador (ou sua combinação linear) aqui, apenas a suavidade da função é a forma mais óbvia de informação a priori que podemos usar.

Amostragem

Falamos sobre funções, mas qualquer dispositivo real não pode medir não apenas um continuum, mas até um conjunto de pontos contáveis. Sempre fazemos medições em um conjunto finito de pontos, portanto somos forçados a executar procedimentos de discretização e transição da equação integral para a matriz. No método de regularização estatística, procedemos da seguinte forma: decomporemos$$$\varphi(x)$ sobre algum sistema de funções $$$\{T_n\}$:

$$

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

Assim, os coeficientes dessa expansão formam algum vetor $$$\vec{\varphi}$que é um vetor no espaço funcional.

Como espaço funcional, podemos pegar o espaço de Hilbert ou, por exemplo, o espaço de polinômios. Além disso, a escolha da base nesses espaços é limitada apenas pela sua imaginação (tentamos trabalhar com as séries trigonométricas de Fourier, poligandra e splines cúbicos).

Então os elementos da matriz$$$K$ calculado como:

$$

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

Onde $$$y_m$- pontos em que as medições foram feitas. Elementos da matriz$$$\Omega$ calcularemos pela fórmula:

$$

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

Onde $$$a$e $$$b$- os limites do intervalo em que a função está definida $$$\varphi(x)$.

Para recalcular erros, use a fórmula de dispersão de uma combinação linear de variáveis ​​aleatórias:

$$

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

Deve-se ter em mente que, em alguns casos, a representação de uma função usando um vetor de dimensão finita leva a uma perda parcial ou alteração de informações. De fato, podemos considerar a algebraização um tipo de regularização, por mais fraca que seja e insuficiente para transformar uma tarefa incorreta em correta. De qualquer forma, agora deixamos de pesquisar$$$\varphi(x)$ para pesquisa de vetor $$$\vec{\varphi}$e na próxima seção, encontramos.

Caso de ruído gaussiano

O caso em que os erros no experimento são distribuídos de acordo com Gauss é notável
, pois uma solução analítica para o nosso problema pode ser obtida. Como as informações e erros a priori têm uma forma gaussiana, seu produto também possui uma forma gaussiana e, em seguida, a feia integral que escrevemos acima é facilmente obtida. A solução e seu erro serão os seguintes:

$$

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

Onde $$$\Sigma$- matriz de covariância da distribuição gaussiana multidimensional, $$$\alpha^*$- o valor mais provável do parâmetro $$$\alpha$, que é determinado a partir da condição de densidade de probabilidade máxima 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})$$

E se eu não tenho erros gaussianos?

A segunda parte do artigo será dedicada a isso, mas por enquanto vamos delinear a essência do problema.

$$

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

O principal problema é que essa terrível integral, primeiramente multidimensional e, em segundo lugar, em limites infinitos. Além disso, é altamente multidimensional, vetor$$$\vec{\varphi}$ pode facilmente ter dimensão $$$m = 30-50$e os métodos de grade para calcular integrais têm complexidade do tipo $$$O(n^m)$, portanto, inaplicável neste caso. Ao usar integrais multidimensionais, a integração Monte Carlo funciona bem.

Além disso, como nossos limites são infinitos, devemos usar o importante método de amostragem, mas precisamos selecionar uma função para amostragem. Para tornar tudo mais automatizado, você deve usar o Monte Carlo da Cadeia de Markov (MCMC) , que pode adaptar independentemente a função de amostragem à integral calculada. Falaremos sobre a aplicação do MCMC no próximo artigo.

Parte prática

A primeira implementação do método de regularização estatística foi escrita nos anos 70 em Algol e foi usada com sucesso para cálculos em física atmosférica. Apesar de ainda termos as fontes manuscritas do algoritmo, decidimos adicionar um pouco de modernismo e fazer uma implementação em Python, e depois em Julia.

Pitão

Instalação

Instale através de pip:

pip install statreg

ou faça o download do código fonte .

Exemplos

Como exemplo, considere como usar um módulo staregpara recuperar dados para uma matriz e equação integral.

Importamos os pacotes científicos necessários.

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

Nós determinamos o sinal verdadeiro, que iremos restaurar.

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

Defina o kernel e a operação de convolução de funções (Nota: np.convolutionespecificamente para matrizes):

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

Geramos os dados medidos e emitimos ruído usando a distribuição 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);

Resolvemos a equação integral

Importamos a classe solver e auxiliar para discretização:

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

Como base funcional para a discretização, usamos splines cúbicos e, como condição adicional, indicamos que a função assume valores zero nas arestas.

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

Resolvemos a equação:

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

Estamos criando um cronograma de soluções:

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

Resolvemos a equação da matriz

Importamos a classe solver e auxiliar para discretização:

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

Para obter matrizes, usamos nossa base funcional, mas é claro que as matrizes podem ser obtidas de qualquer maneira.

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

Resolvemos a equação da matriz:

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

Crie o gráfico:

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

Como mencionamos, o desenvolvimento adicional da técnica requer integração avançada de Monte Carlo. Poderíamos usar algum tipo de módulo em Python (por exemplo, trabalhamos com PyMC3), mas, entre outras coisas, estamos participando de um projeto conjunto com o Instituto Max Planck em Munique.

Este projeto é chamado de Bayesian Analysis Toolkit . Seu objetivo é criar uma estrutura com ferramentas para métodos de análise bayesiana, incluindo principalmente ferramentas para MCMC. Agora, a equipe está trabalhando na segunda versão da estrutura, escrita em Julia (a primeira em C ++ ruim). Uma das tarefas do nosso grupo é demonstrar os recursos dessa estrutura usando o exemplo de regularização estatística, por isso escrevemos uma implementação em 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)

Desta vez, usamos um núcleo diferente, não daremos um passo de integração, mas uma convolução com um gaussiano, que realmente leva um certo "borrão" aos nossos dados:

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)

Da mesma forma, tomamos a base de splines com extremidades fixas:

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)

Exemplo do mundo real

Como exemplo de análise de dados reais, restauraremos o espectro de espalhamento de elétrons de uma mistura hidrogênio-deutério. No experimento, o espectro integral foi medido (ou seja, o número de elétrons acima de uma determinada energia), e precisamos restaurar o espectro diferencial. Para esses dados, o espectro foi reconstruído inicialmente usando o ajuste, para que tenhamos uma base para verificar a correção do nosso algoritmo.

É assim que o espectro integrado inicial se parece:

E assim - o resultado da restauração: A

análise com ajuste tem três desvantagens principais:

• , .
• , .
• .

A regularização estatística evita todos esses problemas e fornece um resultado independente do modelo com erros de medição. A solução obtida pela regularização está de acordo com a curva de ajuste. Observe os dois picos pequenos em 25 e 30 eV. Sabe-se que um pico a 25 eV é formado durante a dispersão dupla e foi restaurado por um ajuste, uma vez que foi claramente especificado na função de ajuste. Um pico de 30 eV pode ser uma anomalia estatística (os erros são bastante grandes neste momento) ou, possivelmente, indica a presença de espalhamento dissociativo adicional.

Conclusões e anúncio da próxima parte

Falamos sobre uma técnica útil que pode ser adaptada a muitas tarefas de análise de dados (incluindo aprendizado de máquina) e para obter um "ajuste" honesto da resposta - a solução mais racional para a equação diante da incerteza causada por erros de medição. Como um bom bônus, obtemos valores para o erro de decisão. Quem quiser participar do desenvolvimento ou aplicar o método de regularização estatística pode contribuir na forma de código em Python, Julia ou em outra coisa.

Na próxima parte, falaremos sobre:

• Usando o MCMC
• Decomposição de Cholesky
• Como exemplo prático, consideramos o uso da regularização para processar um sinal de um modelo de um detector orbital de prótons e elétrons
  

Referências

• Trabalho original Turchin
• A dissertação sobre a aplicação da regularização na física atmosférica
• Nossa publicação sobre a aplicação da regularização estatística
• Implementação Python
• Implementação em Julia
• Implementação no ALGOL-60 , para amantes da arqueologia
• Kit de ferramentas de análise bayesiana

Postado por Mikhail Zeleny , pesquisador do Laboratório de Física Experimental de Métodos da JetBrains Research .