将高斯拟合到测量的峰值

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

scipy.integrate.quad 并没有像您预期的那样进行卷积。 quad(function, lower_bound, upper_bound)[0] 将 return 函数在边界之间的积分的单个值。

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

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

您可能会发现 lmfit (https://lmfit.github.io/lmfit-py/) useful. It has convenient, high-level tools curve-fitting, and simple model functions like Gaussian are built in. It also has a mechanism for fitting a model that is the sum or product of two function, and can even create a model that consists of the convolution of two functions. For example of this, see the examples described in https://lmfit.github.io/lmfit-py/model.html#composite-models-adding-or-multiplying-models.

Answer 2

我觉得你真的很亲近，虽然我不得不承认我不明白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() 来提供背景偏移。好吧，收集更多的数据点 ;).

将高斯拟合到测量的峰值

Fit a Gaussian to measured peak

python

numpy

gaussian

curve-fitting

scipy