获得图的高斯拟合时出错

Question

我不是经验丰富的程序员，需要将我的数据值拟合到高斯图中。下面的代码来自Gaussian fit for Python。

使用Anaconda，得到这个错误信息：

runfile('D:/Anaconda3/gaussian.py', wdir='D:/Anaconda3')

C:\Users\lion\Anaconda3\lib\sitepackages\scipy\optimize\minpack.py:779:OptimizeWarning: Covariance of the parameters could not be estimated category=OptimizeWarning).

拟合的高斯图也不如预期；水平曲线而不是高斯拟合曲线。

图表图像：

任何帮助将不胜感激！

使用的代码：

import pylab as plb
import matplotlib.pyplot as plt

from scipy.optimize import curve_fit
from scipy import asarray as ar,exp


x = ar(range(399))
y = ar([1, 0, 1, 0, 2, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 2, 1, 0, 2, 1, 1, 0, 1, 0, 1, 0, 0, 2, 0, 2, 
        3, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 2, 0, 3, 3, 0, 1, 2, 1, 1, 2, 3, 1, 
        4, 1, 1, 1, 0, 1, 0, 3, 4, 0, 1, 3, 0, 2, 0, 3, 0, 0, 0, 3, 2, 0, 2, 0, 0, 1, 2, 0, 0, 0, 
        3, 2, 1, 0, 1, 3, 4, 2, 4, 1, 2, 1, 1, 2, 0, 2, 2, 6, 2, 4, 2, 0, 1, 2, 2, 3, 4, 6, 2, 3, 
        2, 4, 1, 4, 9, 6, 6, 4, 5, 6, 3, 7, 8, 9, 7, 8, 8, 4, 10, 10, 12, 9, 18, 18, 16, 14, 13, 13, 
        15, 17, 16, 26, 24, 37, 34, 36, 40, 48, 52, 52, 50, 56, 68, 90, 71, 107, 93, 117, 134, 207, 
        200, 227, 284, 287, 337, 379, 449, 471, 626, 723, 848, 954, 1084, 1296, 1481, 1676, 1898, 2024, 
        2325, 2692, 3110, 3384, 3762, 4215, 4559, 5048, 5655, 6092, 6566, 6936, 7513, 8052, 8414, 9016, 
        9303, 9598, 9775, 10100, 10265, 10651, 10614, 10755, 10439, 10704, 10233, 10086, 9696, 9467, 9156, 
        8525, 8200, 7609, 7156, 6678, 6160, 5638, 5227, 4574, 4265, 3842, 3380, 3029, 2767, 2512, 2018, 1856, 1645, 
        1463, 1253, 1076, 943, 787, 711, 588, 512, 448, 361, 304, 303, 251, 190, 185, 154, 134, 114, 105, 
        86, 88, 83, 79, 50, 60, 49, 28, 33, 37, 28, 31, 22, 
        14, 26, 19, 17, 15, 9, 17, 13, 11, 11, 12, 18, 8, 6, 9, 6, 3, 6, 6, 6, 6, 11, 9, 15, 
        3, 3, 1, 2, 3, 2, 6, 3, 4, 3, 4, 5, 3, 1, 1, 2, 1, 0, 4, 3, 2, 3, 1, 3, 3, 4, 0, 3, 5, 0, 
        3, 1, 2, 0, 2, 2, 1, 0, 5, 1, 3, 0, 0, 3, 0, 1, 3, 0, 1, 0, 2, 4, 0, 0, 0, 0, 0, 1, 0, 2, 
        0, 1, 1, 0, 1, 0, 0, 4, 0, 0, 0, 2, 0, 3, 0, 2, 1, 2, 2, 0, 0, 0, 1, 4, 1, 0, 1, 2, 0, 1, 
        1, 1, 1, 2, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 1])

n = len(x)                          #the number of data
mean = sum(x*y)/n                   #note this correction
sigma = sum(y*(x-mean)**2)/n        #note this correction

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

popt,pcov = curve_fit(gaus,x,y,p0=[1,mean,sigma])

plt.plot(x,y,'b+:',label='data')
plt.plot(x,gaus(x,*popt),'ro:',label='fit')
plt.legend()
plt.title('Fig. 3 - Fit for Time Constant')
plt.xlabel('Time (s)')
plt.ylabel('Intensity (Counts)')

plt.show()

Answer 1

在拟合高斯时，正确计算 μ 和 σ 非常重要。您引用了您的参考资料（谢谢！），但没有遵循它。请导入数学并将这两行替换为您的代码：

mean = sum(x * y) / sum(y)
sigma = math.sqrt(sum(y * (x - mean)**2) / sum(y))

（事实证明，如果平均值是介于 -7200 和 10,000 之间的任何数字，拟合仍然可以很好地收敛，对于小于 805 的正 sigma 也是如此。）

Answer 2

您可能会发现 lmfit (http://lmfit.github.io/lmfit-py/) 对这类问题很有用。它具有用于高斯等常见峰形的内置模型，并简化了许多曲线拟合任务。你的例子看起来像这样（跳过数据）：

import matplotlib.pyplot as plt

from lmfit.models import GaussianModel

gmodel = GaussianModel()
params = gmodel.guess(y, x=x)
result = gmodel.fit(y, params, x=x)

print(result.fit_report())

plt.plot(x, y,'b+:',label='data')
plt.plot(x, result.best_fit, 'ro:', label='fit')
plt.legend()
plt.title('Fig. 3 - Fit for Time Constant')
plt.xlabel('Time (s)')
plt.ylabel('Intensity (Counts)')
plt.show()

这将打印出适合的报告

[[Model]]
    Model(gaussian)
[[Fit Statistics]]
    # function evals   = 23
    # data points      = 399
    # variables        = 3
    chi-square         = 1130482.300
    reduced chi-square = 2854.753
    Akaike info crit   = 3177.728
    Bayesian info crit = 3189.695
[[Variables]]
    sigma:       13.0637138 +/- 0.019425 (0.15%) (init= 14.5)
    center:      210.606022 +/- 0.019425 (0.01%) (init= 210.5)
    amplitude:   3.4580e+05 +/- 445.3134 (0.13%) (init= 467842.5)
    fwhm:        30.7626945 +/- 0.045744 (0.15%)  == '2.3548200*sigma'
    height:      10560.0217 +/- 13.59907 (0.13%)  == '0.3989423*amplitude/max(1.e-15, sigma)'
[[Correlations]] (unreported correlations are <  0.100)
    C(sigma, amplitude)          =  0.577

并生成一个显示良好拟合的图。

请注意，高斯的定义略有不同，您的a 将对应于上面列出的height 的值。列为amplitude 的值将是集成区域。