且构网

分享程序员开发的那些事...
且构网 - 分享程序员编程开发的那些事

使用LMFIT将多峰函数拟合到数据集

更新时间:2022-12-23 17:06:00

您应该能够使用内置模型和前缀,如邮件列表中最近有一个关于非常相似主题的讨论. .

You should be able to make use of the built-in models and using prefixes as described in the manual. In addition, there was a recent discussion about a very similar topic on the mailinglist.

您可以执行以下操作.它还不能很好地适应最后一个峰值,但是您可能可以在起始值等方面进行一些调整.而且,由于您的基线并不完全平坦,使用LinearModel而不是ConstantModel可能会有所改善,但是我没有尝试过.

You can do something as shown below. It doesn't fit the last peak very well yet, but you can probably fiddle around a bit with the starting values and such. Moreover, since your baseline isn't completely flat it might improve when you use a LinearModel insteadd of a ConstantModel, but I haven't tried.

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

x, y = np.loadtxt('Peaks.txt', unpack=True)

peaks_in_interval = np.array([43, 159, 191, 296, 435, 544])
number_of_peaks = len(peaks_in_interval)
amplitude = y[peaks_in_interval] / 5
width = np.zeros(number_of_peaks) + 0.1
center = x[peaks_in_interval]

def make_model(num):
    pref = "f{0}_".format(num)
    model = LorentzianModel(prefix = pref)
    model.set_param_hint(pref+'amplitude', value=amplitude[num], min=0, max=5*amplitude[num])
    model.set_param_hint(pref+'center', value=center[num], min=center[num]-0.5, max=center[num]+0.5)
    model.set_param_hint(pref+'sigma', value=width[num], min=0, max=2)
    return model

mod = None
for i in range(len(peaks_in_interval)):
    this_mod = make_model(i)
    if mod is None:
        mod = this_mod
    else:
        mod = mod + this_mod

offset = ConstantModel()
offset.set_param_hint('c', value=np.average(y[-75:]))
mod = mod + offset

out=mod.fit(y, x=x, method='nelder')
plt.interactive(True)
print(out.fit_report())
plt.plot(x, y)
plt.plot(x, out.best_fit, label='best fit')
plt.plot(x, out.init_fit, 'r--', label='fit with initial values')
plt.show()