用lmfit库拟合python中的多峰曲线
提问于 2021-05-19 21:56:34
我有一个数据文件，其中第一列是x-值，第二列是y-值，第三列是y-错误。我想要符合这些数据。我正在学习这里的例子，我的代码是-
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()这段代码给了我下面的情节(不合适)-
我期待着这样的事情-
- 有人能帮我把情节中的数据拟合一下吗？
- 在示例值中，手动定义中心、西格玛和振幅的最小、最大值。有任何方法从数据文件中获取/计算这些值吗？
更新
我试着在评论中使用@mikuszefski建议的峰形。但是，正如图像显示的那样，它也能捕捉到所有的小峰值(噪声)。
是否有办法只为较大的山峰选择价值？
回答 1
Stack Overflow用户
发布于 2021-05-20 09:14:40
你真的必须(“必须”、“被要求”、“在任何情况下绝对”)为所有参数提供合理、有限和合理的初始值。
当你说
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。
您提供的数据和图表显示，您感兴趣的两个峰值发生在大约1.4和1.5的x值上。
所以，中心的值在1左右，你告诉它，这个值在10^5左右，不能低于75,000。在这些限制下，这是最合适的做法。该程序工作正常，没有错误或问题，而且您完全得到了您所要求的。
同样，对于非线性最小二乘问题和曲线拟合，初值总是问题.在任何情况下，它们都是无关紧要的。
也就是说，正如mikuszefski所建议的那样，使用峰值查找算法是一个很好的选择。
旁白:边界应该主要用于约束逻辑/物理。例如，可以合理地说振幅应该是正的。一般的高斯人(或者你的数据)没有内在的要求，要求74999的质心值是荒谬的。所以，不要从这样的界限开始。无限制地开始，尽可能简单地开始。只有在需要时才增加这样的复杂性。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/67614706
复制相关文章
腾讯云开发者
