我试图得到直方图的高斯拟合的拟合误差。我在下面的代码中使用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=Gaussian(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.是否正确地计算了fit中的错误,如果不正确,那么计算它们的正确方法是什么?2.我需要适当的误差来计算一个简化的x-平方值。还有另一种方法,我可以不需要知道在拟合的每个点上的误差,就可以计算齐平方吗?
提前谢谢你!
发布于 2016-02-29 04:09:48
您的主要错误是计算错误值。数组error是对函数Gaussian中系数a,b,c的标准偏差的渐近估计。然后不能输入一个值x和不确定度+/-a,+/-b,+/- c,得到有意义的结果,因为误差大约是a,b&c的平均值,即高斯(x,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)https://stackoverflow.com/questions/35690937
复制相似问题