使用 scipy.curve_fit() 拟合高斯直方图曲线的误差

Question

我正在尝试获取直方图的高斯拟合的拟合误差。我在以下代码中使用 scipy.optimize.curve_fit：

import matplotlib.pylab as plt 
from pylab import exp
import numpy as np 
from scipy import optimize
from math import sqrt

# Fit functions
def Gaussian(x,a,b,c):
    return a * exp(-(x - b)**2.0 / (2 * c**2))

# Generate data from random Guassian distribution
npix = 10200
nbins = int(sqrt(npix))
data = np.random.standard_normal(npix)
print('\n Length of y: %s' % len(data))
n,bins,patches = plt.hist(data,bins=nbins)

# Generate data from bins as a set of points 
bin_size = abs(bins[1]-bins[0])
x =np.linspace(start=bins[0]+bin_size/2.0,stop=bins[-2]+bin_size/2.0,\
num=nbins,endpoint=True)
print min(x),max(x),len(x), np.mean(x)
y = n
y[y==0]= 1e-8


popt, pcov = optimize.curve_fit(Gaussian,x,y) 

# Curve-fit error method
error = [] 
for i in range(len(popt)):
    try:
      error.append( np.absolute(pcov[i][i])**0.5)
    except:
      error.append( 0.00 )
pfit_curvefit = popt
perr_curvefit = np.array(error) 
print('\n Curve-fit Curve fit: %s' % pfit_curvefit)
print('\n Curve-fit Fit errors: %s' % perr_curvefit)

# Plot the fit
x_fit = np.linspace(x[0], x[-1], nbins)
y_gauss = Gaussian(x_fit, *popt)
# y_boot = Gaussian(x_fit, *pfit_bootstrap)
yerr=Gaussian(x_fit,*perr_curvefit)

plt.plot(x_fit, y_gauss,linestyle='--',linewidth=2,\
color='red',label='Gaussian')



plt.xlabel('Pixel Values')
plt.ylabel('Frequency')
plt.title('Npix = %s, Nbins = %s'% (npix,nbins))
plt.legend()
plt.show()

如您所见，我可以 Python 充分拟合直方图数据，没问题。当我尝试计算拟合误差时出现问题

yerr=高斯分布(x_fit,*perr_curvefit)

这似乎是正确的做法，但是当我查看此错误列表时，它看起来毫无意义：

...
0.0
0.0
0.0
0.0
0.0
2.60905702842e-265
2.27384038589e-155
1.02313435685e-74
2.37684931814e-23
0.285080112094
1.76534048255e-08
5.64399121475e-45
9.31623567809e-111
7.93945868459e-206
0.0
0.0
0.0
0.0
0.0
...

我的问题是： 1. 拟合中的误差是否计算正确，如果不正确，计算它们的正确方法是什么。 2. 我需要拟合误差来计算减少的卡方值。有没有另一种方法可以计算卡方，而不必知道拟合中每个点的误差？

提前致谢！

Answer 1

你的主要错误是错误值的计算。您的数组 error 是函数 Gaussian 中系数 a,b,c 的标准差的渐近估计。然后，您不能输入值 x 和不确定性 +/-a、+/-b、+/-c 并获得有意义的结果，因为误差与 a、b 和 c 的平均值有关，即 Gaussian(x, a+/- delta a 等）

如果您不习惯使用 optimize.curve_fit() 并且会使用 optimize.leastsq()，您想要的信息很容易获得。
看到这个 question

替换

popt, pcov = optimize.curve_fit(Gaussian,x,y)

和

def residual(p, x, y):
    return Gaussian(x, *p) - y

initGuess = [1,1,1] # or whatever you want the search to start at
popt, pcov, infodict, mesg, ier = optimize.least_squares(residual,initGuess, args=[x,y], full_output=True)

然后按照解决方案中的说明找到减少的卡方

s_sq = (infodict['fvec']**2).sum()/ (N-n)

使用 scipy.curve_fit() 拟合高斯直方图曲线的误差

Errors on a Gaussian histogram curve fit using scipy.curve_fit()

python

numpy

optimization

scipy

curve-fitting