0

我正在尝试使用 Python 库lmfit拟合数据。在数据中，有两个高斯函数，第二个（LO）高于第一个（TO）。代码如下：

TOmod = GaussianModel(prefix="TO_")
LOmod = GaussianModel(prefix="LO_")
TOmod.set_param_hint('TO_center', min = 265, max = 270)
TOmod.set_param_hint('TO_amplitude', min = 3e3, max = 7e3)
TOmod.set_param_hint('TO_sigma', min = 0.1, max = 5)
LOmod.set_param_hint('LO_center', min = 280, max = 310)
LOmod.set_param_hint('LO_amplitude', min = 2e4, max = 4e4)
LOmod.set_param_hint('LO_sigma', min = 0.1, max = 8)
pars = TOmod.guess(y, x=x)
pars += LOmod.guess(y, x=x)
mod = TOmod + LOmod
out  = mod.fit(y, pars, x=x)
print(out.fit_report())

plt.figure("Fitting")
plt.plot(x, y, "bo")
plt.plot(x, out.init_fit, 'k--')
plt.plot(x, out.best_fit, 'r-')
plt.xlabel(r"k [cm$^{-1}$]")
plt.ylabel(r"Intensity")
plt.xlim(250,325)
plt.ylim(-1e3,3e4)
plt.show("Fitting")

我可以制作这个图表：

lmfit还生成以下报告：

[[Model]]
    (Model(gaussian, prefix='TO_') + Model(gaussian, prefix='LO_'))
[[Fit Statistics]]
    # function evals   = 30
    # data points      = 576
    # variables        = 6
    chi-square         = 2534929941.545
    reduced chi-square = 4447245.511
    Akaike info crit   = 8823.259
    Bayesian info crit = 8849.395
[[Variables]]
    TO_sigma:       0.14748432 +/- 0        (0.00%) (init= 0.1)
    TO_center:      270        +/- 0        (0.00%) (init= 270)
    TO_amplitude:   3000       +/- 0        (0.00%) (init= 3000)
    TO_fwhm:        0.34729903 +/- 0        (0.00%)  == '2.3548200*TO_sigma'
    TO_height:      8114.94309 +/- 0        (0.00%)  == '0.3989423*TO_amplitude/max(1.e-15, TO_sigma)'
    LO_sigma:       0.25064604 +/- 0        (0.00%) (init= 0.1)
    LO_center:      292.364593 +/- 0        (0.00%) (init= 292.3646)
    LO_amplitude:   24876.9938 +/- 0        (0.00%) (init= 20000)
    LO_fwhm:        0.59022631 +/- 0        (0.00%)  == '2.3548200*LO_sigma'
    LO_height:      39595.6185 +/- 0        (0.00%)  == '0.3989423*LO_amplitude/max(1.e-15, LO_sigma)'
[[Correlations]] (unreported correlations are <  0.100)

起初，最终和初始 LO_amplitudes 与报告中的不对应（您可以看到，初始值为 2e4，但在图中至少为 5e4）。怎么来的？

其次，我期待合身会更好，而所有界限都已设定。

编辑 1 - 添加数据

我已经发布了完整的脚本，只有缺少的是加载数据：

x, y = np.genfromtxt("data.txt", unpack=True)

以下是数据：

3.235230000000000246e+02 8.074899999999997817e+02
3.217950000000000159e+02 7.387500000000000000e+02
3.200659999999999741e+02 8.103400000000001455e+02
3.183369999999999891e+02 9.050399999999999636e+02
3.166080000000000041e+02 1.176100000000000364e+03
3.148790000000000191e+02 1.483189999999999600e+03
3.131490000000000009e+02 1.729449999999999818e+03
3.114189999999999827e+02 2.281949999999999818e+03
3.096890000000000214e+02 2.486050000000000182e+03
3.079580000000000268e+02 2.867739999999999782e+03
3.062269999999999754e+02 3.205949999999999818e+03
3.044950000000000045e+02 4.065239999999999782e+03
3.027629999999999768e+02 5.081539999999999964e+03
3.010310000000000059e+02 6.767100000000000364e+03
2.992989999999999782e+02 9.268700000000000728e+03
2.975659999999999741e+02 1.334320000000000073e+04
2.958319999999999936e+02 1.946429999999999927e+04
2.940989999999999895e+02 2.552240000000000146e+04
2.923650000000000091e+02 2.720209999999999854e+04
2.906309999999999718e+02 2.314690000000000146e+04
2.888960000000000150e+02 1.642840000000000146e+04
2.871610000000000014e+02 1.048929999999999927e+04
2.854259999999999877e+02 6.923909999999999854e+03
2.836899999999999977e+02 4.836300000000000182e+03
2.819540000000000077e+02 3.501239999999999782e+03
2.802169999999999845e+02 2.686470000000000255e+03
2.784809999999999945e+02 2.227590000000000146e+03
2.767440000000000282e+02 1.781840000000000146e+03
2.750059999999999718e+02 1.582319999999999709e+03
2.732679999999999723e+02 1.520989999999999782e+03
2.715299999999999727e+02 2.011989999999999782e+03
2.697919999999999732e+02 3.021930000000000291e+03
2.680529999999999973e+02 4.754640000000000327e+03
2.663140000000000214e+02 5.088010000000000218e+03
2.645740000000000123e+02 3.515579999999999927e+03
2.628340000000000032e+02 2.159310000000000400e+03
2.610939999999999941e+02 1.190029999999999745e+03
2.593530000000000086e+02 7.985300000000002001e+02
2.576120000000000232e+02 5.780700000000001637e+02
2.558710000000000093e+02 4.897800000000002001e+02
2.541289999999999907e+02 3.914800000000000182e+02
2.523870000000000005e+02 3.046999999999998181e+02
2.506450000000000102e+02 3.270000000000000000e+02
4

2 回答 2

0

通过遵循我给出的建议并提供初始值而不是施加约束，我得到了相当不错的结果。你试过吗？

LOFWIW，对于更宽的峰值，我使用洛伦兹函数获得了更好的结果。尝试这个：

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

x, y = np.genfromtxt("data.txt", unpack=True)

TOmod = GaussianModel(prefix="TO_")
LOmod = LorentzianModel(prefix="LO_")
mod = TOmod + LOmod

pars = mod.make_params(TO_center=268, TO_sigma=2., TO_amplitude=50000,
                       LO_center=292, LO_sigma=4., LO_amplitude=50000)

out  = mod.fit(y, pars, x=x)
print(out.fit_report())

对于图表报告之间的所谓差异：

正如您可能在文档中读到的那样（http://lmfit.github.io/lmfit-py/builtin_models.html， 是单位归一化amplitude高斯函数前面的乘数，它不是该函数将采用的最大值。 OTOH, 参数和被报告并将是函数的最大值。类似地，方差不是全宽半最大值，但也被报告。LO_heightTO_heightsigma

于 2017-04-16T18:21:40.710 回答
0

看起来它对我来说应该更好。我会推荐两件事：

  1. 通过设置参数值给出明确的初始猜测（特别是对于TO_center和TO_amplitude那些远离边界的）。

  2. 放宽一些界限-这种适合应该很容易，而且我认为界限并不是真正必要的。

如果这些不能解决问题，请发布数据和完整的脚本。

于 2017-04-16T02:00:36.470 回答