使用NumPy / SciPy进行包括所有误差的线性拟合

Question

使用NumPy / SciPy进行包括所有误差的线性拟合

9

我有很多带有y误差的x-y数据点，需要拟合非线性函数。这些函数在某些情况下可以是线性的，但更通常是指数衰减、高斯曲线等等。SciPy支持使用scipy.optimize.curve_fit进行此类拟合，并且我还可以指定每个点的权重。这给了我加权非线性拟合，非常好用。从结果中，我可以提取参数及其各自的误差。

只有一个小问题：误差仅用作权重，但未包含在误差中。如果我将所有数据点的误差加倍，我会期望结果的不确定性也会增加。因此，我构建了一个测试案例（源代码）来测试这一点。

使用scipy.optimize.curve_fit进行拟合，结果如下：

Parameters: [ 1.99900756  2.99695535]
Errors:     [ 0.00424833  0.00943236]

同样的，但是使用2 * y_err：

Parameters: [ 1.99900756  2.99695535]
Errors:     [ 0.00424833  0.00943236]

同样但是y_err乘以2：

所以你可以看到这些值是相同的。这告诉我算法没有考虑到它们，但我认为这些值应该不同。

我也在这里读到了另一种拟合方法，所以我也尝试使用scipy.odr进行拟合：

Beta: [ 2.00538124  2.95000413]
Beta Std Error: [ 0.00652719  0.03870884]

同样的，但是使用 20 * y_err：

Beta: [ 2.00517894  2.9489472 ]
Beta Std Error: [ 0.00642428  0.03647149]

这些数值略有不同，但我认为这并不影响误差的增加。我认为这只是四舍五入误差或者权重略微不同而已。

有没有一款软件包可以让我拟合数据并获得实际误差？我的书里有相关公式，但如果不必要我就不想自己实现。

我现在了解到另一个问题中的 linfit.py 很好地处理了我的需求。它支持两种模式，第一种正是我需要的。

Fit with linfit:
Parameters: [ 2.02600849  2.91759066]
Errors:     [ 0.00772283  0.04449971]

Same but with 20 * y_err:
Parameters: [ 2.02600849  2.91759066]
Errors:     [ 0.15445662  0.88999413]

Fit with linfit(relsigma=True):
Parameters: [ 2.02600849  2.91759066]
Errors:     [ 0.00622595  0.03587451]

Same but with 20 * y_err:
Parameters: [ 2.02600849  2.91759066]
Errors:     [ 0.00622595  0.03587451]

我应该回答我的问题还是现在关闭/删除它？

- Martin Ueding

也许statsmodels可以做到这一点；我不确定它是否处理一般的曲线拟合。 - Fred Foo

4

不要扔掉你所写的一切——回答它，说不定有人知道更好的方法来完成它。 - eickenberg

一定要用你找到的内容来回答问题（同时感谢你在之前一个关于scipy.odr的回答中的评论）。 - Ffisegydd

在scipy 0.14中，curve_fit有一个名为absolute_sigma的选项，请参考http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.curve_fit.html。在此添加之前，对其含义进行了长时间的讨论。 - Josef

@user333700 我的系统上安装了0.13版本，我认为这个选项还没有加入。这是在0.14版本中才添加的吗？那么我可能需要等待或手动安装它。 - Martin Ueding

3个回答

4

请注意，根据curvefit的文档说明：

sigma：None或N长度序列如果不是None，则此向量将在最小二乘问题中用作相对权重。

关键点在于作为相对权重，因此在第53行的yerr和第57行的2*yerr中，您应该会得到类似的结果，如果不是完全相同的结果。

当您增加实际残差误差时，您会看到协方差矩阵中的值变大。例如，如果我们将函数generate_data()中的y += random更改为y += 5 * random，则会发生这种情况。

Fit with scipy.optimize.curve_fit:
('Parameters:', array([ 1.92810458,  3.97843448]))
('Errors:    ', array([ 0.09617346,  0.64127574]))

与原始结果相比较：

Fit with scipy.optimize.curve_fit:
('Parameters:', array([ 2.00760386,  2.97817514]))
('Errors:    ', array([ 0.00782591,  0.02983339]))

另外请注意，参数估计现在距离(2,3)更远了。这是因为残差误差增加和参数估计的置信区间变大所预期的。

- CT Zhu

你是如何增加实际残差误差的？你使用的是0.14还是scipy或其他什么东西？ - Martin Ueding

没有，我还没有更新到0.14.x版本。我只是在你提供的代码中的generate_data()函数中将y += random更改为y += 5*random。 - CT Zhu

2

简短回答

对于包含y（以及odr情况下的x）不确定性的绝对值：

在scipy.odr情况下，使用输出中odr给出的协方差矩阵的cov = numpy.sqrt(numpy.diag(cov))。
在scipy.optimize.curve_fit情况下，使用absolute_sigma=True标志。

对于相对值（不包括不确定性）：

在scipy.odr情况下，使用输出中的sd值。
在scipy.optimize.curve_fit情况下，使用absolute_sigma=False标志。
像这样使用numpy.polyfit：

p, cov = numpy.polyfit(x, y, 1,cov = True) errorbars = numpy.sqrt(numpy.diag(cov))

详细回答

所有函数中都存在一些未记录的行为。我的猜测是这些函数混合了相对和绝对值。最后，该答案是基于如何处理输出而给出您想要的代码（或不提供？是否有错误？）。此外，curve_fit可能最近已获得了“absolute_sigma”标志？

我的观点在于输出。似乎odr计算标准偏差时没有不确定性，类似于polyfit，但是如果从协方差矩阵中计算标准偏差，则存在不确定性。curve_fit使用absolute_sigma=True标志执行此操作。以下是包含的输出

协方差矩阵cov（0,0）和
cov（1,1）的对角线元素，
从斜率和输出的标准偏差的错误方法，以及
常数的错误方法和
从斜率和输出的标准偏差的正确方法，以及
常数的正确方法

odr: 1.739631e-06 0.02302262 [ 0.00014863 0.0170987 ] [ 0.00131895 0.15173207] curve_fit: 2.209469e-08 0.00029239 [ 0.00014864 0.01709943] [ 0.0004899 0.05635713] polyfit: 2.232016e-08 0.00029537 [ 0.0001494 0.01718643]

请注意，odr和polyfit的标准差完全相同。Polyfit不会将不确定性作为输入，因此当计算标准差时，odr不使用不确定性。协方差矩阵使用它们，如果在odr情况下标准偏差是从协方差矩阵中计算出来的，则存在不确定性，并且如果增加不确定性，则会发生变化。在下面的代码中调整dy将显示出来。

我在这里写这个主要是因为在找出误差限时这一点很重要（scipy所引用的fortran odrpack指南中有一些误导性信息：标准差应该像指南所说的那样是协方差矩阵的平方根，但实际上并不是）。

import scipy.odr
import scipy.optimize
import numpy

x = numpy.arange(200)
y = x + 0.4*numpy.random.random(x.shape)
dy = 0.4

def stddev(cov): return numpy.sqrt(numpy.diag(cov))

def f(B, x): return B[0]*x + B[1]

linear = scipy.odr.Model(f) 
mydata = scipy.odr.RealData(x, y,  sy = dy)
myodr = scipy.odr.ODR(mydata, linear, beta0 = [1.0, 1.0], sstol = 1e-20, job=00000)
myoutput = myodr.run()
cov = myoutput.cov_beta
sd  = myoutput.sd_beta
p   = myoutput.beta 
print 'odr:        ', cov[0,0], cov[1,1], sd, stddev(cov)

p2, cov2 = scipy.optimize.curve_fit(lambda x, a, b:a*x+b, 
                                    x, y, [1,1],
                                    sigma = dy,
                                    absolute_sigma = False,
                                    xtol = 1e-20)

p3, cov3 = scipy.optimize.curve_fit(lambda x, a, b:a*x+b, 
                                    x, y, [1,1],
                                    sigma = dy,
                                    absolute_sigma = True,
                                    xtol = 1e-20)

print 'curve_fit:  ', cov2[0,0], cov2[1,1], stddev(cov2), stddev(cov3)

p, cov4 = numpy.polyfit(x, y, 1,cov = True)
print 'polyfit:    ', cov4[0,0], cov4[1,1], stddev(cov4)

- Juha

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Martin Ueding · Accepted Answer

一种行之有效且能够得到更好结果的方法是引入自助法（bootstrap method）。当给定具有误差的数据点时，采用参数自助法，让每个 x 和 y 值描述一个高斯分布。然后从这些分布中各抽取一个点，并获得一个新的自助样本。进行简单的非加权拟合可得到参数的一个值。

该过程重复进行约300到几千次。最终将获得拟合参数的分布，可以计算平均值和标准偏差以获得值和误差。

另一个有趣的事情是，不会获得单个的拟合曲线作为结果，而是很多条。对于每个插值的 x 值，可以再次计算许多值 f(x, param) 的平均值和标准偏差，从而获得一个误差带。

然后，使用各种适合参数再次执行分析的进一步步骤数百次。如上图所示，这也将考虑拟合参数之间的相关性：尽管对数据进行了对称函数的拟合，但误差带是不对称的。这意味着左侧的插值值具有比右侧更大的不确定性。请保留HTML标签。