使用scipy.optimize.curve_fit进行加权拟合

Question

根据文档，参数sigma可用于在拟合中设置数据点的权重。当参数absolute_sigma=True时，它们“描述”1σ误差。

我有一些带人工正态分布噪声的数据，这些噪声是变化的：

n = 200
x = np.linspace(1, 20, n)
x0, A, alpha = 12, 3, 3

def f(x, x0, A, alpha):
    return A * np.exp(-((x-x0)/alpha)**2)

noise_sigma = x/20
noise = np.random.randn(n) * noise_sigma
yexact = f(x, x0, A, alpha)
y = yexact + noise

如果我想使用curve_fit将带噪声的y拟合到f，那么我应该将sigma设置为什么？文档在这里并不是非常具体，但通常我会使用1/noise_sigma**2作为权重：

p0 = 10, 4, 2
popt, pcov = curve_fit(f, x, y, p0)
popt2, pcov2 = curve_fit(f, x, y, p0, sigma=1/noise_sigma**2, absolute_sigma=True)

然而，它似乎并没有显著提高拟合度。

这个选项只是用来通过协方差矩阵更好地解释适应性不确定性的吗？这两者告诉我什么区别？

In [249]: pcov
Out[249]: 
array([[  1.10205238e-02,  -3.91494024e-08,   8.81822412e-08],
       [ -3.91494024e-08,   1.52660426e-02,  -1.05907265e-02],
       [  8.81822412e-08,  -1.05907265e-02,   2.20414887e-02]])

In [250]: pcov2
Out[250]: 
array([[ 0.26584674, -0.01836064, -0.17867193],
       [-0.01836064,  0.27833   , -0.1459469 ],
       [-0.17867193, -0.1459469 ,  0.38659059]])

- xnx

---

当你说它似乎没有改善适合度时，你期望看到什么？ - Croad Langshan

大群的角马在平原上壮观地奔跑。如果这不行，我认为在“带有sigma”的情况下，RMS适配残差会更好，但实际上更糟（0.64与1.07）。 - xnx

LOL，角马。那么未加权算法不是最小化均方根误差吗（回忆一下我做曲线拟合的日子）？这种情况下，加权只会增加它，对吧？你是在告诉它：“不要太担心这些点，更好地拟合这些其他点，即使以牺牲总均方根误差为代价。” - Croad Langshan

注意：R的nls函数可以使用权重，看起来Python中的“sigma”对应于nls的权重的平方根。 - user3637203

Answer 1

至少在 scipy 版本 1.1.0 中，参数 sigma 应该等于每个参数的误差。具体来说，文档 中写道：

一个一维的 sigma 应该包含 ydata 中误差的标准偏差值。在这种情况下，优化函数是 chisq = sum((r / sigma) ** 2)。

对于你的情况，应该是：

curve_fit(f, x, y, p0, sigma=noise_sigma, absolute_sigma=True)

我查看了源代码，确认了当您以这种方式指定sigma时，它会将((f-data)/sigma)**2最小化。

顺便说一句，当您知道误差时，通常希望最小化的就是观察到数据点data给出模型f的可能性：

L(data|x0,A,alpha) = product over i Gaus(data_i, mean=f(x_i,x0,A,alpha), sigma=sigma_i)

如果您取负对数，则变为（除了不依赖于参数的常数因子）：

-log(L) = sum over i (f(x_i,x0,A,alpha)-data_i)**2/(sigma_i**2)

这只是卡方检验。

我编写了一个测试程序来验证curve_fit是否确实返回了正确的值，并且正确指定了sigma：

from __future__ import print_function
import numpy as np
from scipy.optimize import curve_fit, fmin

np.random.seed(0)

def make_chi2(x, data, sigma):
    def chi2(args):
        x0, A, alpha = args
        return np.sum(((f(x,x0,A,alpha)-data)/sigma)**2)
    return chi2

n = 200
x = np.linspace(1, 20, n)
x0, A, alpha = 12, 3, 3

def f(x, x0, A, alpha):
    return A * np.exp(-((x-x0)/alpha)**2)

noise_sigma = x/20
noise = np.random.randn(n) * noise_sigma
yexact = f(x, x0, A, alpha)
y = yexact + noise

p0 = 10, 4, 2

# curve_fit without parameters (sigma is implicitly equal to one)
popt, pcov = curve_fit(f, x, y, p0)
# curve_fit with wrong sigma specified
popt2, pcov2 = curve_fit(f, x, y, p0, sigma=1/noise_sigma**2, absolute_sigma=True)
# curve_fit with correct sigma
popt3, pcov3 = curve_fit(f, x, y, p0, sigma=noise_sigma, absolute_sigma=True)

chi2 = make_chi2(x,y,noise_sigma)

# double checking that we get the correct answer
xopt = fmin(chi2,p0,xtol=1e-10,ftol=1e-10)

print("popt  = %s, chi2 = %.2f" % (popt,chi2(popt)))
print("popt2 = %s, chi2 = %.2f" % (popt2, chi2(popt2)))
print("popt3 = %s, chi2 = %.2f" % (popt3, chi2(popt3)))
print("xopt  = %s, chi2 = %.2f" % (xopt, chi2(xopt)))

输出结果如下：

popt  = [ 11.93617403   3.30528488   2.86314641], chi2 = 200.66
popt2 = [ 11.94169083   3.30372955   2.86207253], chi2 = 200.64
popt3 = [ 11.93128545   3.333727     2.81403324], chi2 = 200.44
xopt  = [ 11.93128603   3.33373094   2.81402741], chi2 = 200.44

从结果可以看出，当您将sigma=sigma作为curve_fit函数的参数指定时，chi2确实被正确地最小化了。

至于为什么改进效果不是“更好”，我不太确定。我的唯一猜测是，如果没有指定sigma值，则隐含地假定它们相等，并且在拟合很重要的数据部分（峰值）上，误差“大约”相等。

回答您的第二个问题，不仅仅是用sigma选项来更改协方差矩阵输出，它实际上改变了正在最小化的内容。

- user545424

---

为什么不把sigma称为噪声？ - Kornpob Bhirombhakdi

如果您知道噪声项，那么您可以将其从数据中减去，然后您就有了一个完美的信号，甚至不需要拟合任何东西。对于真实数据，您通常知道误差的标准偏差，但您不知道每个数据点的实际误差，这就是为什么需要进行拟合的原因。 - user545424