使用 lmfit 进行高斯拟合答案

【问题标题】：Gaussian Fit using lmfit使用 lmfit 进行高斯拟合
【发布时间】：2021-05-14 17:56:09
【问题描述】：

我正在尝试拟合高斯。我尝试使用 OriginPro 和 Python 来适应。 OriginPro 的拟合比通过 Python 获得的要好，我想用 Python 来做。

我使用的代码是：

import numpy as np
import matplotlib.pyplot as plt
import lmfit

data=data
mod=lmfit.models.GaussianModel()
x=np.arange(140,510,1)
y= np.array(np.log(data[140:510]), dtype=np.float32)
idx= (np.isfinite(y))
pars = mod.guess(y[idx], x=x[idx])
out = mod.fit(y[idx], pars, x=x[idx])
plt.plot(x[idx], out.best_fit, 'r', label='Gaussian Model')
plt.plot(x,np.log(data[140:510]))
plt.show()

数据：

[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., 34., 33., 30.,
       25., 47., 36., 37., 33., 41., 42., 40., 49., 47., 44., 35.,  0.,
        0.,  0., 30., 32., 37., 36., 50., 40., 30., 41., 29., 33., 34.,
       42., 41., 39., 43., 34., 30., 27., 48., 38., 20., 53., 38., 32.,
       41., 40., 45., 39., 34., 33., 41., 32., 41., 41., 30., 42., 41.,
       33., 62., 36., 47., 44., 35., 44., 36., 40., 35., 43., 38., 36.,
       48., 50., 31., 46., 50., 43., 50., 40., 35., 39., 50., 47., 46.,
       35., 39., 33., 49., 43., 56., 45., 56.,  0.,  0.,  0.,  0., 42.,
       41., 43., 46., 40., 40., 38., 40., 39., 49., 33., 39., 57., 51.,
       51., 43., 47., 55., 67., 68., 44., 48.,  0.,  0.,  0.,  0.,  0.,
       48., 66., 55., 60., 63., 58., 72., 53., 62., 55., 63., 64., 60.,
       61., 56., 58., 71., 51., 74., 58., 69., 61., 50., 59., 66., 57.,
       56., 68., 63., 53., 70., 73., 49., 65., 56., 62., 45., 55.,  0.,
        0.,  0.,  0., 66., 53., 52., 52., 61., 42., 53., 49., 81., 67.,
       63., 69., 58., 80., 74., 51., 55., 65., 68., 56., 71., 58., 59.,
       84., 54., 49., 46., 69., 57., 54., 74., 68., 60., 56., 51., 58.,
       47., 57.,  0.,  0.,  0.,  0., 65., 47., 59., 39., 34., 57., 45.,
       43., 57., 61., 39., 52., 48., 42., 53., 41.,  0.,  0.,  0.,  0.,
       48., 35., 53., 44., 46., 52., 50., 51., 44., 55., 42., 58., 53.,
       42.,  0.,  0.,  0.,  0., 41., 39., 47., 46., 46., 46., 41., 44.,
       49., 48., 47., 43., 46., 38., 47., 39., 29., 53., 40., 45., 43.,
       41., 23., 41., 39., 46., 37., 42., 24., 40., 46., 39., 37., 49.,
       47., 37., 28., 36., 34., 34., 43., 41., 48., 42., 44., 45., 43.,
       45., 36., 33., 48., 39., 54., 35., 30., 42., 34., 51., 39., 50.,
       36., 34., 50., 37., 48., 40., 39., 21., 41., 40., 42., 44., 48.,
       42., 42., 39., 43., 57., 37., 36., 38., 34., 50., 40., 30., 29.,
       38., 44., 43., 47., 35., 37., 44., 44., 32., 36., 41., 44., 28.,
       43., 36., 37., 34., 37.,  0.,  0.,  0.,  0., 54., 34., 45.,  0.,
        0.,  0.,  0.,  0., 44., 32., 28., 34., 42., 35., 40., 45., 42.,
       47., 39., 35., 52., 34., 55., 38.,  1.]

fit_report 的输出：

In [37]:out.fit_report()
Out[37]: "[[Model]]\n    Model(gaussian)\n[[Fit Statistics]]\n    # fitting method   = leastsq\n    # function evals   = 17\n    # data points      = 333\n    # variables        = 3\n    chi-square         = 12.5923742\n    reduced chi-square = 0.03815871\n    Akaike info crit   = -1084.59201\n    Bayesian info crit = -1073.16758\n[[Variables]]\n    amplitude:  3590.97370 +/- 125.195985 (3.49%) (init = 781.4034)\n    center:     317.738344 +/- 3.57625327 (1.13%) (init = 322.4771)\n    sigma:      361.440601 +/- 13.6760525 (3.78%) (init = 181.5)\n    fwhm:       851.127556 +/- 32.2046420 (3.78%) == '2.3548200*sigma'\n    height:     3.96355944 +/- 0.01657428 (0.42%) == '0.3989423*amplitude/max(1e-15, sigma)'\n[[Correlations]] (unreported correlations are < 0.100)\n    C(amplitude, sigma) =  0.997"

和我使用的猜测函数有关系吗？

【问题讨论】：

不要将代码作为图片发布
见math.stackexchange.com/q/1292889/295281。我在一个名为 scikit-guess 的包中实现了这个方法，你可以在 pypi 上找到它。它仍然是实验性的
您的代码不够完整，无法评论可能出了什么问题。提供一个完整的、最小的示例并显示结果（来自 fit_report()）以及绘图。通过快速查看图像，您的 OriginPro 代码似乎包含了一个常量偏移量以及一个高斯，因此您可能也希望将其包含在 Python 代码中。
@MNewville 现在我已经添加了所有相关的代码，你想让我也显示 2theta 中的数据吗？我还添加了 fit_report() 并且情节已经存在。
@Pranav45 这不可能是所有代码。您忽略了导入，尽管这些可能很明显。更重要的是，theta_up 是未定义的。结果，没有人能诊断出你的代码做了什么。见stackoverflow.com/help/minimal-reproducible-example。而且：您的模型不包括常量（可能是变化的？）偏移量。您的数据不需要吗？

标签： python curve-fitting lmfit

【解决方案1】：

在 cmets 中给出的提示似乎不够。您需要在模型中包含偏移量。高斯函数远离峰值强度为 0 - 我不知道 OriginPro 在做什么，但很明显，它不仅仅是一个高斯函数。

尝试制作一个高斯 + 常数的模型，如下所示：

import numpy as np
import matplotlib.pyplot as plt
import lmfit

data = np.array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              34., 33., 30., 25., 47., 36., 37., 33., 41., 42., 40.,
              49., 47., 44., 35., 0., 0., 0., 30., 32., 37., 36., 50.,
              40., 30., 41., 29., 33., 34., 42., 41., 39., 43., 34.,
              30., 27., 48., 38., 20., 53., 38., 32., 41., 40., 45.,
              39., 34., 33., 41., 32., 41., 41., 30., 42., 41., 33.,
              62., 36., 47., 44., 35., 44., 36., 40., 35., 43., 38.,
              36., 48., 50., 31., 46., 50., 43., 50., 40., 35., 39.,
              50., 47., 46., 35., 39., 33., 49., 43., 56., 45., 56.,
              0., 0., 0., 0., 42., 41., 43., 46., 40., 40., 38., 40.,
              39., 49., 33., 39., 57., 51., 51., 43., 47., 55., 67.,
              68., 44., 48., 0., 0., 0., 0., 0., 48., 66., 55., 60.,
              63., 58., 72., 53., 62., 55., 63., 64., 60., 61., 56.,
              58., 71., 51., 74., 58., 69., 61., 50., 59., 66., 57.,
              56., 68., 63., 53., 70., 73., 49., 65., 56., 62., 45.,
              55., 0., 0., 0., 0., 66., 53., 52., 52., 61., 42., 53.,
              49., 81., 67., 63., 69., 58., 80., 74., 51., 55., 65.,
              68., 56., 71., 58., 59., 84., 54., 49., 46., 69., 57.,
              54., 74., 68., 60., 56., 51., 58., 47., 57., 0., 0., 0.,
              0., 65., 47., 59., 39., 34., 57., 45., 43., 57., 61.,
              39., 52., 48., 42., 53., 41., 0., 0., 0., 0., 48., 35.,
              53., 44., 46., 52., 50., 51., 44., 55., 42., 58., 53.,
              42., 0., 0., 0., 0., 41., 39., 47., 46., 46., 46., 41.,
              44., 49., 48., 47., 43., 46., 38., 47., 39., 29., 53.,
              40., 45., 43., 41., 23., 41., 39., 46., 37., 42., 24.,
              40., 46., 39., 37., 49., 47., 37., 28., 36., 34., 34.,
              43., 41., 48., 42., 44., 45., 43., 45., 36., 33., 48.,
              39., 54., 35., 30., 42., 34., 51., 39., 50., 36., 34.,
              50., 37., 48., 40., 39., 21., 41., 40., 42., 44., 48.,
              42., 42., 39., 43., 57., 37., 36., 38., 34., 50., 40.,
              30., 29., 38., 44., 43., 47., 35., 37., 44., 44., 32.,
              36., 41., 44., 28., 43., 36., 37., 34., 37., 0., 0., 0.,
              0., 54., 34., 45., 0., 0., 0., 0., 0., 44., 32., 28.,
              34., 42., 35., 40., 45., 42., 47., 39., 35., 52., 34.,
              55., 38., 1.])

data = data[140:510]
ipos = np.where(data>0)[0]
data = np.log(data[ipos])
x    = np.arange(140,510,1.0)[ipos]

mod  = lmfit.models.GaussianModel()  + lmfit.models.ConstantModel()

pars = mod.make_params(c=data.mean(),
                       center=x.mean(),
                       sigma=x.std(),
                       amplitude=x.std()*data.ptp())

out = mod.fit(data, pars, x=x)

print(out.fit_report())

plt.plot(x, out.best_fit, 'r', label='Gaussian Model')
plt.plot(x, data)
plt.show()

这将打印一份报告

[Model]]
    (Model(gaussian) + Model(constant))
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 48
    # data points      = 333
    # variables        = 4
    chi-square         = 9.50683177
    reduced chi-square = 0.02889615
    Akaike info crit   = -1176.19189
    Bayesian info crit = -1160.95932
[[Variables]]
    amplitude:  58.7086539 +/- 5.27950257 (8.99%) (init = 153.376)
    center:     302.462920 +/- 2.79518589 (0.92%) (init = 322.6336)
    sigma:      48.1997878 +/- 3.66967747 (7.61%) (init = 106.8759)
    c:          3.64131210 +/- 0.01679060 (0.46%) (init = 3.796494)
    fwhm:       113.501824 +/- 8.64142989 (7.61%) == '2.3548200*sigma'
    height:     0.48592258 +/- 0.02641186 (5.44%) == '0.3989423*amplitude/max(1e-15, sigma)'
[[Correlations]] (unreported correlations are < 0.100)
    C(amplitude, c)     = -0.832
    C(amplitude, sigma) =  0.798
    C(sigma, c)         = -0.664

还有一个情节

【讨论】：

有什么方法可以访问变量的不确定性吗？例如，幅度的不确定性是指 5.27950257
最佳拟合参数及其不确定性将在out.params中，这是一个以参数名称为键的字典，值是具有value、stderr、vary等属性的Parameter对象， min、max、correl 等。所以使用out.params['amplitude'].value、out.params['amplitude'].stderr 等