python - curve_fit 和 scipy.odr 的比较 - 绝对西格玛

Question

我目前想要拟合 x 和 y 中有错误的数据，我使用 scipy.odr 包来获得我的结果。我只是想知道如何正确使用错误 sx 和 sy。所以这是一个例子。

假设我针对不同的电流 I 测量电压 V。 所以我测量了 1V，误差为 +- 0.1V 等等。因此，如果我假设电流测量没有错误，我可以按如下方式使用 spipy.curve_fit。 Absolute_sigma 设置为 True，因为我的绝对误差为 +- 0.1V。

所以我得到:

from scipy.optimize import curve_fit
y=[1,2.5,3,4]
x=[1,2,3,4]
yerr=[0.1,0.1,0.1,0.1]

def func(x, a, b):
    return a*x+b

popt, pcov = curve_fit(func, x, y,sigma=yerr,absolute_sigma=True)

print(popt)
print(np.sqrt(np.diag(pcov)))

[0.95 0.25]
[0.04472136 0.12247449]

在第二步中，我想使用电流和电压均有误差的 odr 包。 根据文档，它应该按如下方式使用:sx 和 sy 是我的测量数据的误差。 因此，如果我对 sx 使用非常小的误差，我应该假设得到与 curve_fit 相似的结果。

from scipy.odr import *
x_err = [0.0000000001]*x.__len__()
y_err = yerr

def linear_func(p, x):
   m, c = p
   return m*x+c

linear_model = Model(linear_func)

data = RealData(x, y, sx=x_err, sy=y_err)

odr = ODR(data, linear_model, beta0=[0.4, 0.4])

out = odr.run()

out.pprint()

Beta: [0.94999996 0.24999994]
Beta Std Error: [0.13228754 0.36228459]

但是如您所见，结果与上面 curve_fit absolute_sigma=True 的结果不同。

将相同的数据与 curve_fit 和 absolute_sigma=False 一起使用会导致与 ODR-Fit 相同的结果。

popt, pcov = curve_fit(func, x, y,sigma=yerr,absolute_sigma=False)

print(popt)
print(np.sqrt(np.diag(pcov)))

[0.95 0.25]
[0.13228756 0.36228442]

所以我猜 ODR-Fit 并没有像 curve_fit 和 absolute_sigma=True 那样真正处理我的绝对误差。有什么办法可以做到这一点，还是我遗漏了什么？

Answer 1

curve_fit() 中的选项 absolute_sigma=True 给出了 实协方差矩阵 在这个意义上 np.sqrt( np.diag(pcov)) 产生 1-sigma 标准偏差作为定义的误差，例如 Numerical Recipes，即它可以有一个单位，如 meter 左右。 kmpfit 附带一个非常有用的摘要

但是，ODR 给出了从缩放协方差矩阵导出的标准误差。看这里，How to compute standard error from ODR results?或下面的示例。

通过 ODR 执行这种缩放，使得减少的 chi2 - 使用以相同方式缩放的输入权重计算 - 产生大约 1。或来自 curve_fit docs :

This constant is set by demanding that the reduced chisq for the optimal parameters popt when using the scaled sigma equals unity.

现在剩下的问题是:sd_beta 到底是什么意思？ 这是标准错误，参见例如Standard error of mean versus standard deviation (可能存在两者大小相同的情况。请参阅下面的评论) 另见前面的讨论:https://mail.python.org/pipermail/scipy-user/2013-February/034196.html

现在，通过使用缩减的 chi2 缩放 pcov，对于参数误差，获得与 a) curve_fit(..., absolute_sigma=False) 和 b 相同的输出) ODR a priory 相对错误:

# calculate chi2
residuals = (y - func(x, *popt))
chi_arr =  residuals / yerr
chi2_red = (chi_arr**2).sum() / (len(x)-len(popt))
print('red. chi2:\t\t\t\t', chi2_red)
print('np.sqrt(np.diag(pcov) * chi2_red):\t', np.sqrt(np.diag(pcov) * chi2_red))

产量:

red. chi2:                           8.75
np.sqrt(np.diag(pcov) * chi2_red):   [0.13228756 0.36228442]

请注意，但是，对于 curve_fit(..., absolute_sigma=True) 和 ODR，协方差矩阵 pcov 来自 curve_fit 和ODR 输出的 cov_beta 仍然相同。在这里，来回重新缩放的缺点变得显而易见!

相对误差 现在当然意味着，如果缩放输入误差，则相对误差的大小保持不变:

# scaled uncertainty
yerr = np.asfarray([0.1, 0.1, 0.1, 0.1]) * 2

popt, pcov = curve_fit(func, x, y, sigma=yerr, absolute_sigma=False)

print('popt:\t\t\t', popt)
print('np.sqrt(np.diag(pcov)):\t', np.sqrt(np.diag(pcov)))

与相同的输出:

popt:                    [0.95 0.25]
np.sqrt(np.diag(pcov)):  [0.13228756 0.36228442]