【发布时间】：2016-10-12 20:42:49
【问题描述】：
我正在尝试使用 scipy 的 curve_fit 拟合分布。我试图拟合一个单分量指数函数，结果几乎是一条直线（见图）。我还尝试了两个分量的指数拟合，它似乎工作得很好。两个分量只是意味着方程的一部分以不同的输入参数重复。无论如何，这里是单分量拟合函数：
def Exponential(Z,w0,z0,Z0):
z = Z - Z0
termB = (newsigma**2 + z*z0) / (numpy.sqrt(2.0)*newsigma*z0)
termA = (newsigma**2 - z*z0) / (numpy.sqrt(2.0)*newsigma*z0)
return w0/2.0 * numpy.exp(-(z**2 / (2.0*newsigma**2))) * (numpy.exp(termA**2)*erfc(termA) + numpy.exp(termB**2)*erfc(termB))
装修完毕
fitexp = curve_fit(Exponential,newx,y2)
然后我尝试了一些东西，只是为了尝试一下。我取了两个分量拟合的两个参数，但在计算中没有用到。
def ExponentialNew(Z,w0,z0,w1,z1,Z0):
z = Z - Z0
termB = (newsigma**2 + z*z0) / (numpy.sqrt(2.0)*newsigma*z0)
termA = (newsigma**2 - z*z0) / (numpy.sqrt(2.0)*newsigma*z0)
return w0/2.0 * numpy.exp(-(z**2 / (2.0*newsigma**2))) * (numpy.exp(termA**2)*erfc(termA) + numpy.exp(termB**2)*erfc(termB))
突然之间就可以了。
现在，我的条件是。为什么？如您所见，拟合的计算绝对没有区别。它只是获得两个未使用的额外变量。这不应该得到相同的结果吗？
@Andras_Deak 一个实际的例子：
from scipy.special import erfc
import numpy
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
#setup data
x = [-58.,-54.,-50.,-46.,-42.,-38.,-34.,-30.,-26.,-22.,-18.,-14.,-10.,-6.,-2.,2.,6.,10.,14.,18.,22.,26.,30.,34.,38.,42.,46.,50.,54.,58.]
y = [23.06763817, 16.89802085, 17.83258379, 16.63446237, 13.81878965, 12.97965839, 14.30451789, 16.98288216, 22.26811491, 28.56756908, 33.06990344, 38.59842098, 54.19860393, 86.37381604, 137.47253315, 199.49724512, 238.66047662, 219.89405445, 160.68820199, 103.88901303, 65.92405727, 43.84596266, 31.5395342, 25.9610156, 22.71683709, 18.06740651, 13.85362374, 11.12867065, 10.36502799, 11.31855619]
y_err = [17.9823065, 4.13684885, 1.66490726, 2.4109372, 2.93359141, 1.9701747, 3.19214881, 3.65593012, 2.89089074, 3.58922121, 4.25505348, 4.72728874, 6.77736567, 11.3888196, 21.87771722, 39.0087495, 56.6910311, 51.7592369, 26.39750958, 10.62678862, 7.85893395, 8.11741621, 7.91731416, 7.07739132, 5.41818744, 6.11286843, 8.27070757, 7.85323065, 4.26885499, 0.9047867]
#function to fit
def Exponential2(Z, w0, z0, w1, z1, Z0):
z = Z - Z0
s = 3.98098937586
a = z**2 / (2.0*s**2)
b = (s**2 + z*z0) / (numpy.sqrt(2.0)*s*z0)
c = (s**2 - z*z0) / (numpy.sqrt(2.0)*s*z0)
d = (s**2 + z*z1) / (numpy.sqrt(2.0)*s*z1)
e = (s**2 - z*z1) / (numpy.sqrt(2.0)*s*z1)
return w0/2.0 * numpy.exp(-a) * (numpy.exp(c**2)*erfc(c) + numpy.exp(b**2)*erfc(b)) + w1/2.0 * numpy.exp(-a) * (numpy.exp(e**2)*erfc(e) + numpy.exp(d**2)*erfc(d))
#derive and set initial guess
ymaxpos = x[numpy.where(y==numpy.max(y))[0]]
p0_2 = [numpy.max(y),5,numpy.max(y)/2.0,20,ymaxpos]
#fit
fitexp2 = curve_fit(Exponential2,x,y,p0=p0_2,sigma=y_err)
#get results
w0err = numpy.sqrt(numpy.diag(fitexp2[1]))[0]
z0err = numpy.sqrt(numpy.diag(fitexp2[1]))[1]
w1err = numpy.sqrt(numpy.diag(fitexp2[1]))[2]
z1err = numpy.sqrt(numpy.diag(fitexp2[1]))[3]
w0 = fitexp2[0][0]
z0 = fitexp2[0][1]
w1 = fitexp2[0][2]
z1 = fitexp2[0][3]
Z0 = fitexp2[0][4]
#new x array for smoother curve
smoothx = numpy.arange(-58,59,0.1)
y2 = Exponential2(smoothx,w0,z0,w1,z1,Z0)
print 'Exponential 2: w0: '+str(w0.round(3))+' +/- '+str(w0err.round(3))+' \t z0: '+str(z0.round(3))+' +/- '+str(z0err.round(3))+' \t w1: '+str(w1.round(3))+' +/- '+str(w1err.round(3))+' \t\t z1: '+str(z1.round(3))+' +/- '+str(z1err.round(3))
#plot
fig = plt.figure()
ax = fig.add_subplot(111)
ax.errorbar(x,y,y_err,fmt='o',markersize=2,label='data')
ax.plot(smoothx,y2,label='fit',color='red')
ax.grid()
ax.legend()
plt.show()
如您所见，情节看起来不错，但返回的值 z1 完全不现实。
Exponential 2: w0: 312.608 +/- 36.764 z0: 8.263 +/- 1.158 w1: 12.689 +/- 9.138 z1: 1862257.883 +/- 45201809883.8
【问题讨论】：
-
这是一个我不熟悉的相当复杂的模型。它不对应两个驼峰吗？您的输入数据在左侧有一个巨大的驼峰和一个很小的驼峰，但应该很难找到后者。我怀疑你的模型中的一个驼峰远离你感兴趣的数据，所以基本上你正在为你的数据拟合一个驼峰（因此是无意义的参数）。您可以尝试拟合这个大驼峰，从数据中减去它，然后再将另一个拟合到剩余部分。
-
实际上，减去该主峰会留下一个噪声数据集，该数据集似乎没有明显值得注意的特征（在误差线内）。您确定该模型非常适合您的问题吗？可能是部分信息不相关/不必要（例如将两个高斯的总和拟合到一个峰值）。
-
是的。这是沿星系短轴的光分布，应该由该模型描述。然而，我注意到，如果我去掉前 5 个数据点，结果会更可信。也许curve_fit 无法处理这第一个碰撞......
-
这很有可能。如果该凹凸明显独立于大凹凸，则您确实应该适当地忽略它。问题是最大的不确定性远大于左侧凸块上的预期噪声，即使该凸块在误差范围内高于基线。