对它们的不正确逆问题的统计正则化。Turchin(第1部分)

哈Ha!今天我们要告诉您的是JetBrains Research的一部分,核物理实验方法实验室的工作内容 你问JetBrains在哪里,核物理在哪里。我们是基于对Kotlin的热爱而达成的,尽管在这篇文章中我们不会谈论他。我们的团队专注于为科学家开发数据分析,建模和编写软件,因此致力于与IT公司的合作和知识共享。 在本文中,我们要讨论由V.F. Turchin在XX世纪70年代提出的普及的统计正则化方法,以及以Python和Julia的代码形式实现的方法。





演示将非常详细,因此那些了解反问题的人可以直接阅读示例并阅读本文中的理论


问题的发生:为什么任何人都应该正规化?


如果足以忽略不计,那么实验中的任何测量都可以描述如下:有某种设备可以捕获过程的频谱或信号,并根据测量结果显示一些数字。作为研究人员,我们的任务是查看这些数字并了解设备结构,以了解被测频谱或信号是什么。也就是说,面对所谓的逆问题如果您从数学上进行想象,我们将得到以下方程式(顺便说一句,称为第一种Fredholm方程式):

f(y)=abdxK(x,y)φ(x)

实际上,此等式描述了以下内容:我们的测量设备在这里用其硬件功能表示 K(x,y)对研究的频谱或其他输入信号起作用 φ使研究人员观察输出信号 f(y)研究人员的目标是恢复信号φ 由著名 f(y)K(x,y)您也可以用矩阵形式表示此表达式,用向量和矩阵替换函数:

fm=Kmnφn

信号重建似乎不是一项艰巨的任务,因为Fredholm方程和线性方程组(甚至超定)都具有精确的解决方案。因此,让我们尝试一下。让被测信号描述为两个高斯之和:

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

作为一种工具,我们采用最简单的积分器-使用Heaviside函数将信号转换为累积和的矩阵:

Kmn=θ(xmyn)

图中显示了被测信号的类型和我们的设备以及测量结果。

重要的是,任何实际测量都存在误差,因此我们将通过添加正常噪声来稍微破坏我们的结果,从而给出5%的测量误差。

我们将通过最小二乘法恢复信号:

φ=(KTK)1KTf

结果我们得到:

实际上,我们可以完成这篇文章,再次使自己确信面对现实的残酷而无情的数学理想主义方法是无奈的,然后去吸烟。

但是首先,让我们找出导致这种失败发生在我们身上的原因是什么?显然,关键是测量误差,但是它们会产生什么影响?事实是,雅克·哈达玛(Jacques Hadamard)(向柯西-哈达玛公式加了破折号的那个人)将所有任务分为正确摆放的任务和不正确摆放的任务

记住经典:“寻找解决方案(如果有的话)是没有意义的。它是关于如何处理没有解决方案的任务。这是一个深层的根本问题……”-我们不会谈论正确的任务,而立即处理不正确的任务。幸运的是,我们已经遇到了这样的问题:上面写的Fredholm方程是一个不正确的逆问题- 即使输入数据有无限小的波动(甚至我们的测量误差也远非无穷小),以精确的分析方式获得的方程的解也可能与真实的解任意不同

您可以在A.N.院士的经典著作的第一章中阅读此陈述的证据。 Tikhonova“解决不适定问题的方法。”本书提供了处理不正确任务的技巧,但其中概述的技术在Turchin方法中已解决了许多缺陷。但是首先,我们描述处理不正确任务的一般原则:遇到此类任务该怎么办?

由于任务本身无法为我们提供任何东西,因此有必要犯下小罪行:用数据补充任务,使其变得正确,换句话说,输入一些关于任务的先验信息*(此过程称为任务正则化)。与基于参数化正则化函数引入的经典Tikhonov方法相反,Turchin方法则使用贝叶斯方法进行调用。

统计正则化的理论描述


战略


我们根据数理统计表述我们的问题:根据一个著名的实现 f(我们在实验中进行了测量),我们需要评估参数的值 φ功能性S^ 计算 φ基于 f我们将称为策略为了确定哪种策略更理想,我们引入了二次损失函数实际损失函数可以是任何函数,为什么我们选择二次函数?因为任何接近最小值的损失函数都可以用二次函数来近似:

L(φ,S^[f])=||φ^S^[f])||L2,

哪里 φ^-最佳解决方案。然后,我们选择的策略的损失由风险函数确定

RS^[f](φ)E[L(φ,S^[f])]=L(φ,S^[f])P(f|φ)df.

这里 P(f|φ)确定集合的概率密度,在该集合上进行平均损失。该合奏由假设的多次测量重复形成。f 给定的 φ通过这种方式,P(f|φ) -这是我们已知的概率密度 f在实验中获得。

根据贝叶斯方法,建议考虑φ作为具有先验概率密度随机变量 P(φ)表达我们解决问题各种解决方案可靠性P(φ)根据实验之前存在的信息确定。然后,基于最小化后验风险来选择最佳策略

rS^(φ)EφEf[L(φ,S^[f])|φ]

在这种情况下,最佳策略是众所周知的:

S^[f]=E[φ|f]=φP(φ|f)dφ,

后密度在哪里 P(φ|f)由贝叶斯定理确定:

P(φ|f)=P(φ)P(f|φ)dφP(φ)P(f|φ)

这种方法将使我们能够确定所得解决方案的方差(相关函数):

D(x1,x2)=E[φ(x1)S^[f](x1)][φ(x2)S^[f](x2)]

因此,我们通过引入先验密度获得了解决问题的最佳方案 P(φ)我们能谈谈功能世界吗φ(x)由先验密度给出的?

如果这个问题的答案是否定的,那么我们将不得不接受所有可能的情况。φ(x)同样可能,并返回不规则解决方案。因此,我们必须肯定地回答这个问题。

这正是统计正则化方法的本质所在-通过引入有关以下信息的先验信息来对解决方案进行正则化φ(x)如果研究人员已经具有任何先验信息(先验密度P(φ)),他只需要计算积分并获得答案即可。

如果没有此类信息,下一段将描述研究人员可能拥有的最少信息,以及如何使用它来获取正规化的解决方案。

先验信息


正如英国科学家所表明的那样,他们喜欢在世界其他地区与众不同。此外,如果数学家会问到有关此操作合法性的问题,物理学家乐观地认为,自然法则是通过“良好”函数(即平滑函数)来描述的。

换句话说,它使它更平滑φ(x)先验概率密度更高。因此,让我们尝试基于平滑度引入先验概率。为此,我们记得,先验信息的引入是对世界的某种暴力,迫使自然法则以对我们方便的方式看待。

应将这种暴力行为减至最少,并通过引入先验概率密度,有必要在信息φ(x)包含在 P(φ)是最小的。形式化上述过程,我们基于函数的平滑度导出先验密度的形式。为此,我们将寻找一个有条件的信息极值:

I[P(φ)]=lnP(φ)P(φ)dφmin

在以下条件下:

  1. 光滑度条件 φ(x)Ω是表征函数平滑度的某个矩阵。然后我们要求达到一定的平滑度功能值:

    (φ,Ωφ)P(φ)dφ=ω

    细心的读者应该问一个有关确定参数值的问题。
    ω答案将在本文中进一步给出。
  2. 单位概率归一化: P(φ)dφ=1
    在这种情况下,以下功能将提供最少的功能:

    Pα(φ)=αRg(Ω)/2detΩ1/2(2π)N/2exp(12(φ,αΩφ))

    参数 α连接 ω但是,由于我们没有有关平滑功能特定值的信息,因此弄清楚它的连接方式毫无意义。那该怎么办α, 你问。这里向您展示了三种方式:

    1. 匹配参数值 α手动操作,从而实际进行Tikhonov的正则化
    2. 平均(积分)所有可能 α假设所有可能 α同样可能
    3. 选择最可能的 α通过其后验概率密度 P(α|f)这种方法是正确的,因为如果实验数据包含有关以下内容的足够信息,则它可以很好地近似积分α

第一种情况对我们来说并不重要。在第二种情况下,我们必须在这里计算这样一个丑陋的积分:

φi=dφφiP(f|φ)dαP(α)αRg(Ω)2exp(α2(φ,Ωφ))dφP(f|φ)dαP(α)αRg(Ω)2exp(α2(φ,Ωφ))

对于第三种情况,我们可以在实验中解析获得高斯噪声的积分值(将在本节中进行讨论)。

还应该指出,我们没有在任何地方使用过Ω是平滑度运算符。实际上,在这里我们可以使用任何其他运算符(或它们的线性组合),只是函数的平滑度是我们可以使用的先验信息的最明显形式。

采样


我们讨论了功能,但是任何真实的设备都不能测量一个连续体,甚至不能测量一组点。我们总是在一组有限的点上进行测量,因此我们必须执行从积分方程到矩阵方程的离散化和转换过程。在统计正则化方法中,我们进行如下操作:我们将分解φ(x) 在某些功能系统上 {Tn}

φ(x)=nφnTn(x).

因此,这个扩展的系数形成了一个向量 φ这是函数空间中的向量。

作为函数空间,我们可以采用希尔伯特空间,例如多项式空间。此外,在这些空间中选择基础仅受您的想象力限制(我们尝试使用三角傅里叶级数,波利甘德拉和三次样条)。

然后矩阵的元素K 计算为:

Kmn=(K^Tn(x))(ym),

哪里 ym-进行测量的点。矩阵元素Ω 我们将通过以下公式计算:

Ωij=ab(dpTi(x)dx)(dpTj(x)dx)dx,

哪里 ab-定义函数的时间间隔的边界 φ(x)

要重新计算误差,请使用随机变量线性组合的色散公式:

D[φ(x)]=D[nφnTn(x)]=i,jφiφjcov(Ti(x),Tj(x)).

应该记住,在某些情况下,使用有限维向量表示函数会导致部分信息丢失或更改。实际上,我们可以将代数视为一种正则化,但是它微不足道且不足以将不正确的任务变成正确的任务。但是,无论如何,我们现在已经不再使用搜索φ(x) 矢量搜索 φ在下一节中,我们将找到它。

高斯噪声案例


根据高斯分布实验中的误差的情况非常显着
,因为可以获得对我们问题的解析解决方案。由于先验信息和错误具有高斯形式,因此它们的乘积也具有高斯形式,因此我们上面编写的丑陋积分很容易采用。解决方案及其错误如下:

φ=(KTΣ1K+αΩ)1KTΣ1Tf

Σφ=(KTΣ1K+αΩ)1,

哪里 Σ-多维高斯分布的协方差矩阵, α-参数的最可能值 α,这取决于最大后验概率密度的条件:

P(α|f)=CαRg(Ω)2|(KTΣ1K+αΩ)1|exp(12fTΣ1KT(KTΣ1K+αΩ)1KTΣ1Tf)


如果我没有高斯误差?


本文的第二部分将专门讨论这一点,但现在让我们概述问题的实质。

φi=dφφiP(f|φ)dαP(α)αRg(Ω)2exp(α2(φ,Ωφ))dφP(f|φ)dαP(α)αRg(Ω)2exp(α2(φ,Ωφ))


主要问题是这种可怕的积分,首先是多维的,其次是无限的。而且,它是高度多维的矢量φ 可以轻松拥有尺寸 m=3050,并且用于计算积分的网格方法具有类型复杂度 O(nm)因此,在这种情况下不适用。当采用多维积分时,蒙特卡洛积分效果很好。

此外,由于我们的限制是无限的,因此必须使用重要的采样方法,但是随后我们需要选择一个函数进行采样。为了使一切自动化,您应该使用马尔可夫链蒙特卡洛(MCMC),它可以独立地使采样函数适应所计算的积分。我们将在下一篇文章中讨论MCMC的应用。

实践部分


统计正则化方法的第一个实现是在70年代在Algol上写的,并成功地用于大气物理学中的计算。尽管我们仍然拥有该算法的手写源,但我们还是决定添加一些现代性,并在Python和Julia上实现。

蟒蛇


安装

通过安装pip

pip install statreg

或下载源代码

例子

例如,考虑如何使用模块stareg为矩阵和积分方程恢复数据。

我们导入必要的科学软件包。

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

我们确定真实信号,并将其恢复。

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


定义内核和函数的卷积运算(注意:np.convolution专门针对数组):

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

我们生成测量数据并使用正态分布对其进行噪声处理:

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



我们求解积分方程

我们导入求解器和辅助类以进行离散化:

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

作为离散化的功能基础,我们使用三次样条曲线,作为附加条件,我们指示该函数的边缘取零值。

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

我们求解方程:

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

我们正在制定解决方案时间表:

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



我们求解矩阵方程

我们导入求解器和辅助类以进行离散化:

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

要获得矩阵,我们使用我们的功能基础,但是很显然,可以以任何方式获得矩阵。

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

我们求解矩阵方程:

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

建立图表:

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



朱莉亚


如前所述,该技术的进一步发展需要先进的蒙特卡洛集成。我们可以在Python中使用某些模块(例如,我们与PyMC3一起使用),但是除其他外,我们正在与慕尼黑的马克斯·普朗克研究所共同参与一个项目。

该项目称为贝叶斯分析工具包。其目标是创建一个包含用于贝叶斯分析方法的工具的框架,主要包括用于MCMC的工具。现在,团队正在研究该框架的第二个版本,该版本是用Julia编写的(第一个版本是用不良C ++编写的)。我们小组的任务之一是使用统计正则化示例来演示此框架的功能,因此我们在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)


这次我们使用不同的核心,我们将不执行积分步骤,而是与高斯进行卷积运算,这实际上会使我们的数据产生一定的“模糊”:

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)


同样,我们以固定端的样条曲线为基础:

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)



真实的例子

作为实际数据分析的一个示例,我们将还原氢-氘混合物的电子散射光谱。在实验中,对积分光谱进行了测量(即,电子数高于一定能量),我们需要恢复微分光谱。对于此数据,首先使用拟合来重建光谱,以便我们有一个基础来检查算法的正确性。

这就是初始积分谱的外观:

等等-恢复的结果:

拟合分析具有三个主要缺点:

  • , .
  • , .
  • .

统计正则化避免了所有这些问题,并提供了与模型无关的结果,并带有测量误差。通过正则化得到的解与拟合曲线非常吻合。注意在25和30 eV处的两个小峰。众所周知,在两次散射过程中会在25 eV处形成一个峰,并且由于在拟合函数中明确指定了峰,因此通过拟合将其恢复。30 eV的峰值可能是统计异常(此时的误差非常大),或者可能表明存在其他解离散射。

结论和下一部分的公告


我们向您介绍了一种有用的技术,该技术可以适应许多数据分析任务(包括机器学习),并诚实地“解决”问题-面对由测量误差引起的不确定性,方程的最合理解决方案。作为一个不错的奖励,我们获得决策错误的值。那些想要参与开发或应用统计正则化方法的人可以以Python,Julia或其他形式的代码形式做出贡献。

在下一部分中,我们将讨论:

  • 使用MCMC
  • 胆固醇分解
  • 作为一个实际例子,我们考虑使用正则化处理来自质子和电子轨道探测器模型的信号

参考文献



发贴者米哈伊尔沉降值,研究员在核物理实验方法的实验室JetBrains的研究

All Articles