在 python 中使用 lmfit 库进行多峰曲线拟合

Question

我有一个datafile，其中第一列是 x 值，第二列是 y 值，第三列是 y 误差。我想拟合数据。我正在遵循here 的示例，我的代码是-

import matplotlib.pyplot as plt
import numpy as np

from lmfit.models import ExponentialModel, GaussianModel


file='sample-data.txt'
dat = np.loadtxt(file)
x = dat[:, 0]
y = dat[:, 1]



exp_mod = ExponentialModel(prefix='exp_')
pars = exp_mod.guess(y, x=x)

gauss1 = GaussianModel(prefix='g1_')
pars.update(gauss1.make_params())

pars['g1_center'].set(value=105000, min=75000, max=125000)
pars['g1_sigma'].set(value=150000, min=30000)
pars['g1_amplitude'].set(value=2000000, min=100000)

gauss2 = GaussianModel(prefix='g2_')
pars.update(gauss2.make_params())

pars['g2_center'].set(value=155000, min=125000, max=175000)
pars['g2_sigma'].set(value=150000, min=30000)
pars['g2_amplitude'].set(value=2000000, min=100000)

mod = gauss1 + gauss2 + exp_mod

init = mod.eval(pars, x=x)
out = mod.fit(y, pars, x=x)

print(out.fit_report(min_correl=0.5))

fig, axes = plt.subplots(1, 2, figsize=(12.8, 4.8))
axes[0].plot(x, y, 'b')
axes[0].plot(x, init, 'k--', label='initial fit')
axes[0].plot(x, out.best_fit, 'r-', label='best fit')
axes[0].legend(loc='best')

comps = out.eval_components(x=x)
axes[1].plot(x, y, 'b')
axes[1].plot(x, comps['g1_'], 'g--', label='Gaussian component 1')
axes[1].plot(x, comps['g2_'], 'm--', label='Gaussian component 2')
axes[1].plot(x, comps['exp_'], 'k--', label='Exponential component')
axes[1].legend(loc='best')

plt.show()

这段代码给了我下面的情节（拟合不起作用）-

我期待这样的事情-

1. 谁能帮我拟合图中的数据？
2. 还在example 值、最小值、最大值中手动定义了中心、西格玛和幅度。有没有办法从数据文件中获取/计算这些值？

更新

我尝试使用@mikuszefski 在评论中建议的find_peaks。但它也拾取了所有小峰值（噪声），如图所示。

有没有办法只为较大的峰选择值？

Answer 1

您确实必须（“必须”、“被要求”、“在所有情况下都绝对”）为所有参数提供合理、有限且合理的初始值。

当你说

pars['g1_center'].set(value=105000, min=75000, max=125000)
pars['g1_sigma'].set(value=150000, min=30000)
pars['g1_amplitude'].set(value=2000000, min=100000)

gauss2 = GaussianModel(prefix='g2_')
pars.update(gauss2.make_params())

pars['g2_center'].set(value=155000, min=125000, max=175000)
pars['g2_sigma'].set(value=150000, min=30000)
pars['g2_amplitude'].set(value=2000000, min=100000)

你（字面意思是“字面意思”）告诉程序高斯 #1 应该从中心值 105000 开始，并且在任何情况下都不能超过 [75000, 125000]。

您提供的数据和图表显示，您感兴趣的两个峰值出现在 x 值约为 1.4 和 1.5 处。

因此，中心的值约为 1，而您告诉它该值约为 10^5，并且不能低于 75,000。这是在这些限制下的最佳选择。该程序运行正常，没有错误或问题，您得到的正是您所要求的。

同样，对于非线性最小二乘问题和曲线拟合，初始值总是很重要。没有任何情况下它们无关紧要。

也就是说，使用 mikuszefski 建议的寻峰算法是一个不错的选择。

旁白：边界应该主要用于约束逻辑/物理。例如，可以合理地说幅度应该是正的。一般来说，高斯（或者可能是您的数据）没有什么内在的要求质心值 74999 是不合理的。所以，不要从这样的界限开始。开始没有界限，尽可能简单。仅在需要时才添加这种复杂性。