如何正确包含用python拟合的不确定性
我试图用python中的y个不确定性来拟合一些数据点。 数据在Python中被标记为x，y和yerr。 我需要在loglog规模上对这些数据进行线性拟合。 作为一个参考，如果适合的结果是正确的，我比较蟒蛇的结果与Scidavis的结果
我尝试了curve_fit
def func(x, a, b):
return np.exp(a* np.log(x)+np.log(b))
popt, pcov = curve_fit(func, x, y,sigma=yerr)
以及kmpfit
def funcL(p, x):
a,b = p
return ( np.exp(a*np.log(x)+np.log(b)) )
def residualsL(p, data):
a,b=p
x, y, errorfit = data
return (y-funcL(p,x)) / errorfit
a0=1
b0=0.1
p0 = [a0,b0]
fitterL = kmpfit.Fitter(residuals=residualsL, data=(x,y,yerr))
fitterL.parinfo = [{}, {}]
fitterL.fit(params0=p0)
当我试图用没有不确定性的数据（即设置yerr = 1）来拟合数据时，一切正常，结果与scidavis中的数据完全相同。 但是如果我对数据文件的不确定性置之不理，我会得到一些令人不安的结果。 在Python中，我得到ie a = 0.86和scidavis a = 0.14。 我读了一些关于错误被包含为权重的内容。 我是否必须改变任何事情，才能正确计算出合适的身材？ 或者我做错了什么？
编辑：这里是一个数据文件的例子（x，y，yerr）
3.942387e-02 1.987800e+00 5.513165e-01
6.623142e-02 7.126161e+00 1.425232e+00
9.348280e-02 1.238530e+01 1.536208e+00
1.353088e-01 1.090471e+01 7.829126e-01
2.028446e-01 1.023087e+01 3.839575e-01
3.058446e-01 8.403626e+00 1.756866e-01
4.584524e-01 7.345275e+00 8.442288e-02
6.879677e-01 6.128521e+00 3.847194e-02
1.032592e+00 5.359025e+00 1.837428e-02
1.549152e+00 5.380514e+00 1.007010e-02
2.323985e+00 6.404229e+00 6.534108e-03
3.355974e+00 9.489101e+00 6.342546e-03
4.384128e+00 1.497998e+01 2.273233e-02
结果是：
in python:
without uncertainties: a=0.06216 +/- 0.00650 ; b=8.53594 +/- 1.13985
with uncertainties: a=0.86051 +/- 0.01640 ; b=3.38081 +/- 0.22667
in scidavis:
without uncertainties: a = 0.06216 +/- 0.08060; b = 8.53594 +/- 1.06763
with uncertainties: a = 0.14154 +/- 0.005731; b = 7.38213 +/- 2.13653
我一定是误解了一些东西。 您发布的数据看起来不像
f(x,a,b) = np.exp(a*np.log(x)+np.log(b))
红线是scipy.optimize.curve_fit的结果，绿线是scidavis的结果。
我的猜测是，两种算法都不会朝着合适的方向收敛，所以结果不匹配也就不足为奇了。
我无法解释scidavis是如何找到它的参数的，但根据我理解的定义， scipy正在寻找具有比scidavis更低的最小二乘残差的参数：
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as optimize
def func(x, a, b):
return np.exp(a* np.log(x)+np.log(b))
def sum_square(residuals):
return (residuals**2).sum()
def residuals(p, x, y, sigma):
return 1.0/sigma*(y - func(x, *p))
data = np.loadtxt('test.dat').reshape((-1,3))
x, y, yerr = np.rollaxis(data, axis = 1)
sigma = yerr
popt, pcov = optimize.curve_fit(func, x, y, sigma = sigma, maxfev = 10000)
print('popt: {p}'.format(p = popt))
scidavis = (0.14154, 7.38213)
print('scidavis: {p}'.format(p = scidavis))
print('''
sum of squares for scipy: {sp}
sum of squares for scidavis: {d}
'''.format(
sp = sum_square(residuals(popt, x = x, y = y, sigma = sigma)),
d = sum_square(residuals(scidavis, x = x, y = y, sigma = sigma))
))
plt.plot(x, y, 'bo', x, func(x,*popt), 'r-', x, func(x, *scidavis), 'g-')
plt.errorbar(x, y, yerr)
plt.show()
产量
popt: [ 0.86051258 3.38081125]
scidavis: (0.14154, 7.38213)
sum of squares for scipy: 53249.9915654
sum of squares for scidavis: 239654.84276
上一篇: How to correctly include uncertainties in fitting with python
下一篇: NAME and WM