线性拟合，包括 NumPy/SciPy 的所有错误答案

【问题标题】：Linear fit including all errors with NumPy/SciPy线性拟合，包括 NumPy/SciPy 的所有错误
【发布时间】：2016-12-04 20:21:54
【问题描述】：

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

只有一个警告：错误仅用作权重，但不包含在错误中。如果我将所有数据点的误差加倍，我预计结果的不确定性也会增加。所以我建立了一个测试用例（source code）来测试这个。

适合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]

相同但有 2 * y_err:

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

我在这里也读到了另一种拟合方法，所以我也尝试使用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]

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

【问题讨论】：

也许 statsmodels 可以做到这一点；我不确定它是否可以处理一般曲线拟合。
不要扔掉你写的所有东西——回答它，谁知道呢，也许有人知道更好的方法。
用你发现的东西肯定回答你的问题（感谢你对我之前讨论scipy.odr的答案之一的评论）。
在 scipy 0.14 docs.scipy.org/doc/scipy-0.14.0/reference/generated/… 中查看 absolute_sigma 选项 curve_fit docs.scipy.org/doc/scipy-0.14.0/reference/generated/… 在添加此内容之前，对此的含义进行了长时间的讨论。
@user333700 我的系统上安装了 0.13，但我认为还没有这个选项。刚刚在 0.14 中添加？然后我可能不得不等待或手动安装它。

标签： python numpy scipy

【解决方案1】：

请注意，来自curvefit的文档：

sigma : 无或 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) 相差甚远，正如我们预期的那样，残差误差增加和参数估计的置信区间更大。

【讨论】：

您是如何增加实际残留误差的？您使用的是 0.14 还是 scipy 或其他的？
不，已经更新到 0.14.x。我刚刚在您链接到的代码中将函数generate_data() 中的y += random 更改为y += 5*random。

【解决方案2】：

一种行之有效且实际效果更好的方法是引导方法。当给出有错误的数据点时，使用参数引导程序并让每个x 和y 值描述一个高斯分布。然后从这些分布中的每一个中抽取一个点并获得一个新的自举样本。执行简单的未加权拟合会为参数提供一个值。

这个过程会重复大约 300 到几千次。最终会得到拟合参数的分布，其中可以采用均值和标准差来获得值和误差。

另一件巧妙的事情是，我们不会因此获得单一的拟合曲线，而是获得很多拟合曲线。对于每个内插的x 值，我们可以再次取许多值f(x, param) 的平均值和标准差，并获得一个误差带：

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

【讨论】：

【解决方案3】：

简答

对于在 y 中包含不确定性的绝对值（对于 odr 情况，在 x 中）：

在scipy.odr 的情况下使用stddev = numpy.sqrt(numpy.diag(cov)) 其中 cov 是 odr 在输出中给出的协方差矩阵。
在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)

【讨论】：