将高斯函数拟合到测量的峰值

Question

我有一个小的谱峰值，并且我正在试图拟合一个高斯函数。我在网上搜索了一个例子，并将代码与我所做的代码混合在一起。

wveleng=[ 639.188 639.454 639.719 639.985 640.25 640.516 640.781 641.046 
     641.312 641.577] 
    counts=[ 778. 1613.8 12977.4 32990. 33165.2 13171. 2067.2 900.8 
     788.8 747.8] 

我的第一个代码如下

def gaus(x,a,mu,sigma): 
    return a*exp(-(x-mu)**2/(2*sigma**2)) 

    a=ydata.max() 
x0=ydata.mean() 
sigm=ydata.std() 

mean = sum(ydata*xdata)/len(ydata) 
sigma = np.sqrt(sum(ydata*(xdata-mean)**2)/len(ydata)) 

#print(ydata.max()) 
popt, pcov = curve_fit(Gauss, xdata,ydata,maxfev=991,p0=[a,x0,sigm])  
#gmodel = Model(Gauss) 
#result = gmodel.fit(ydata, x=xdata, a=ydata.max(),x0=ydata.mean(),sigm=ydata.std()) 
print(popt) 
#plt.scatter(xdata,ydata,label='data points') 
#plt.plot(xdata, result.best_fit, 'r-') 
#popt, pcov = curve_fit(gauss, xdata, ydata,p0=[ydata.max(), ydata.mean(), ydata.std()]) 
xx = np.linspace(639,642, 10) 
plt.plot(xx, gauss(xdata, *popt), 'r-', label='fit') 

与情节，我得到以下。

我认为这与最初的猜测参数

我觉得这更紧凑，更适合更好地为我的第二代码做。

def gauss(x, a, x0, sigma): 
    return a * np.exp(-(x - x0) ** 2/(2 * sigma ** 2)) 

ydata = np.array([778.,1613.8,12977.4,32990.,33165.2,13171.,2067.2,900.8,788.8,747.8]) 

xx = np.arange(639,642, 100) 
xdata=np.array([639.188,639.454,639.719,639.985,640.250,640.516,640.781,641.046,641.312,641.577]) 


#plt.plot(xdata, ydata, 'bo', label='data') 
def Gauss(x, a, x0, sigm): 
    return a * np.exp(-(x - x0)**2/(2 * sigm**2)) 

gmodel = Model(Gauss) 
result = gmodel.fit(ydata, x=xdata, a=ydata.max(),x0=ydata.mean(),sigm=ydata.std()) 
plt.scatter(xdata,ydata,label='data points') 
plt.plot(xdata, result.best_fit, 'r-') 

我得到和第一种方法完全相同的结果。有没有一种方法比数据本身更适合多点

Answer 1

我认为你真的很接近，但我不得不承认我不明白xx是什么意思。您一定希望数据拟合（ydata）和自变量（xdata）的长度相同。

我想你现在正在运行到的主要问题是，你的初始猜测都不是很好，而且你会得到良好的结果

result = gmodel.fit(ydata, x=xdata, a=ydata.max(), x0=xdata.mean(), sigma=xdata.std()) 

（与xdata而不是ydata控制初始值x0和sigma）。

或许更重要的是一些健全检查添加到参数值的范围，与

params = gmodel.make_params(a=ydata.max(), 
          x0=xdata.mean(), 
          sigma=xdata.std()) 
params['x0'].min = min(xdata) 
params['x0'].max = max(xdata) 
params['sigma'].max = 5 
result = gmodel.fit(ydata, params, x=xdata) 

最后，使用内置的模型，如GaussianModel将同时报告sigma和fwhm也amplitude（即高斯的积分）和height（高斯将采用的最大值）。所以，这个脚本：

import numpy as np 
from lmfit.models import GaussianModel 
import matplotlib.pyplot as plt 

ydata = np.array([778.,1613.8,12977.4,32990.,33165.2,13171.,2067.2,900.8,788.8,747.8]) 
xdata = np.array([639.188,639.454,639.719,639.985,640.250,640.516,640.781,641.046,641.312,641.577]) 

gmodel = GaussianModel() 
params = gmodel.make_params(amplitude=ydata.max(), 
          center=xdata.mean(), 
          sigma=xdata.std()) 
result = gmodel.fit(ydata, params, x=xdata) 

print(result.fit_report()) 
plt.scatter(xdata,ydata,label='data points') 
plt.plot(xdata, result.best_fit, 'r-') 
plt.show() 

会打印出

[[Model]] 
    Model(gaussian) 
[[Fit Statistics]] 
    # function evals = 27 
    # data points  = 10 
    # variables  = 3 
    chi-square   = 2360971.771 
    reduced chi-square = 337281.682 
    Akaike info crit = 129.720 
    Bayesian info crit = 130.628 
[[Variables]] 
    sigma:  0.27525232 +/- 0.004505 (1.64%) (init= 0.7623906) 
    center:  640.119396 +/- 0.004490 (0.00%) (init= 640.3828) 
    amplitude: 25633.2702 +/- 362.5571 (1.41%) (init= 33165.2) 
    fwhm:  0.64816968 +/- 0.010608 (1.64%) == '2.3548200*sigma' 
    height:  37152.0777 +/- 525.0717 (1.41%) == '0.3989423*amplitude/max(1.e-15, sigma)' 
[[Correlations]] (unreported correlations are < 0.100) 
    C(sigma, amplitude)   = 0.579 

，并给予相当好看的配合。对于高级练习，我建议尝试添加一个ConstantModel（）以提供背景偏移量。那么，并收集更多的数据点;）。

Answer 2

scipy.integrate.quad没有进行卷积，正如你所期望的那样。 quad(function, lower_bound, upper_bound)[0]将为边界之间的函数积分返回单个值。

OTOH，curve_fit(func, ...)需要一个模型函数值的数组，它抱怨它有一个浮点数，而不是一个ndarray。

也许你打算做curve_fit(vfunc, ...)？

您可能会觉得lmfit（https://lmfit.github.io/lmfit-py/）有用。它具有方便的高级工具曲线拟合功能，并且内置了高斯等简单模型函数。它还具有拟合作为两个函数的总和或乘积的模型的机制，甚至可以创建一个由两个函数的卷积。例如，请参阅https://lmfit.github.io/lmfit-py/model.html#composite-models-adding-or-multiplying-models中描述的示例。