我正在尝试对某些数据拟合经过指数修改的高斯函数。数据位于顶部。
我有以下代码:
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from scipy.special import erfc
bins = [-46.82455,-46.41738,-46.01021,-45.60304,-45.19587,-44.7887,-44.38153,-43.97436,-43.56719,-43.16002,-42.75285,-42.34568,-41.93851,-41.53134,-41.12417,-40.717,-40.30983,-39.90266,-39.49549,-39.08832,-38.68115,-38.27398,-37.86681,-37.45964,-37.05247,-36.6453,-36.23813,-35.83096,-35.42379,-35.01662,-34.60945,-34.20228,-33.79511,-33.38794,-32.98077,-32.5736,-32.16643,-31.75926,-31.35209,-30.94492,-30.53775,-30.13058,-29.72341,-29.31624,-28.90907,-28.5019,-28.09473,-27.68756,-27.28039,-26.87322,-26.46605,-26.05888,-25.65171,-25.24454,-24.83737,-24.4302,-24.02303,-23.61586,-23.20869,-22.80152,-22.39435,-21.98718,-21.58001,-21.17284,-20.76567,-20.3585,-19.95133,-19.54416,-19.13699,-18.72982,-18.32265,-17.91548,-17.50831,-17.10114,-16.69397,-16.2868,-15.87963,-15.47246,-15.06529,-14.65812,-14.25095,-13.84378,-13.43661,-13.02944,-12.62227,-12.2151,-11.80793,-11.40076,-10.99359,-10.58642,-10.17925,-9.77208,-9.36491,-8.95774,-8.55057,-8.1434,-7.73623,-7.32906,-6.92189,-6.51472,-6.10755,-5.70038,-5.29321,-4.88604,-4.47887,-4.0717,-3.66453,-3.25736,-2.85019,-2.44302,-2.03585,-1.62868,-1.22151,-0.81434,-0.40717,0.0,0.40717,0.81434,1.22151,1.62868,2.03585,2.44302,2.85019,3.25736,3.66453,4.0717,4.47887,4.88604,5.29321,5.70038,6.10755,6.51472,6.92189,7.32906,7.73623,8.1434,8.55057,8.95774,9.36491,9.77208,10.17925,10.58642,10.99359,11.40076,11.80793,12.2151,12.62227,13.02944,13.43661,13.84378,14.25095,14.65812,15.06529,15.47246,15.87963,16.2868,16.69397,17.10114,17.50831,17.91548,18.32265,18.72982,19.13699,19.54416,19.95133,20.3585,20.76567,21.17284,21.58001,21.98718,22.39435,22.80152,23.20869,23.61586,24.02303,24.4302,24.83737,25.24454,25.65171,26.05888,26.46605,26.87322,27.28039,27.68756,28.09473,28.5019,28.90907,29.31624,29.72341,30.13058,30.53775,30.94492,31.35209,31.75926,32.16643,32.5736,32.98077,33.38794,33.79511,34.20228,34.60945,35.01662,35.42379,35.83096,36.23813,36.6453,37.05247,37.45964,37.86681,38.27398,38.68115,39.08832,39.49549,39.90266,40.30983,40.717,41.12417,41.53134,41.93851,42.34568,42.75285,43.16002,43.56719,43.97436,44.38153,44.7887,45.19587,45.60304,46.01021,46.41738]
counts = [0.00000000e+00,0.00000000e+00,9.82318271e-04,1.96463654e-03,7.85854617e-03,9.82318271e-03,1.27701375e-02,1.47347741e-02,1.76817289e-02,2.75049116e-02,3.14341847e-02,4.32220039e-02,5.79567780e-02,6.77799607e-02,9.43025540e-02,1.29666012e-01,1.48330059e-01,1.87622790e-01,2.07269155e-01,2.54420432e-01,3.00589391e-01,3.33005894e-01,4.03732809e-01,4.72495088e-01,5.22593320e-01,5.99214145e-01,6.34577603e-01,7.04322200e-01,8.18271120e-01,8.58546169e-01,9.26326130e-01,9.65618861e-01,9.35166994e-01,9.76424361e-01,9.39096267e-01,1.00000000e+00,9.67583497e-01,9.36149312e-01,9.13555992e-01,9.38113949e-01,8.35952849e-01,8.31041257e-01,8.33988212e-01,7.54420432e-01,7.17092338e-01,6.12966601e-01,6.22789784e-01,5.37328094e-01,4.76424361e-01,4.35166994e-01,3.89980354e-01,3.53634578e-01,3.47740668e-01,3.51669941e-01,2.87819253e-01,2.67190570e-01,3.04518664e-01,2.60314342e-01,2.70137525e-01,2.65225933e-01,3.06483301e-01,2.72102161e-01,2.61296660e-01,2.57367387e-01,2.45579568e-01,2.25933202e-01,2.28880157e-01,2.21021611e-01,2.23968566e-01,1.95481336e-01,1.80746562e-01,1.56188605e-01,1.53241650e-01,1.23772102e-01,1.47347741e-01,1.26719057e-01,8.93909627e-02,7.17092338e-02,8.84086444e-02,6.28683694e-02,6.97445972e-02,6.58153242e-02,4.61689587e-02,4.51866405e-02,4.22396857e-02,3.92927308e-02,3.43811395e-02,2.45579568e-02,3.53634578e-02,2.94695481e-02,2.16110020e-02,3.63457760e-02,1.96463654e-02,2.35756385e-02,2.84872299e-02,2.55402750e-02,2.06286837e-02,1.86640472e-02,3.33988212e-02,2.25933202e-02,2.65225933e-02,6.87622790e-03,1.66994106e-02,1.17878193e-02,1.08055010e-02,1.37524558e-02,5.89390963e-03,9.82318271e-03]
def exp_mod_gauss(x,m,s,l):
y = 0.5*l*np.exp(0.5*l*(2*m+l*s*s-2*x))*erfc((m+l*s*s-x)/(np.sqrt(2)*s))
return y
#l=Lambda,s=Sigma,m=Mu
bins=np.asarray(bins,dtype='float')
counts=np.asarray(counts,dtype='float')
popt,pcov = curve_fit(exp_mod_gauss,bins,counts,p0=[-3.5,2.8736,0.1548])
fitted_func = exp_mod_gauss(bins,popt[0],popt[1],popt[2])
#fitted_func = exp_mod_gauss(bins,-3.5,0.1548) #used for manual example
plt.plot(bins,'o',markersize=1) #plot actual counts
plt.plot(bins,fitted_func/max(fitted_func),'-') #plot fitted func/scaled
plt.show()
如果按照编写的方式使用scipy拟合运行代码,则会得到以下结果:
显然不是很好。
但是,如果我注释掉使用curve_fit参数的fit_func行,并使用我在初始猜测中提供的参数(-3.5、2.876、0.1548),则会得到以下结果:
因此,即使我为curve_fit提供最初的猜测,这基本上就是我要寻找的答案,但是它失败了。通过在Matlab中执行完全相同的过程,我得到了很好的拟合参数,但是我不想使用Matlab。我想使用Python。
有人知道这里发生了什么吗?
非常感谢。