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

Question

我正在尝试使用 LMFIT 库进行多洛伦兹拟合，但它不起作用，我什至理解我的语法made 是完全错误的，但我没有任何新的想法。

我的问题是这样的:我有一个很长的光谱，其中包含多组峰，但这些组中峰的数量不是恒定的，因此有时我只有 1 个峰值，但有时可能有 8 个甚至 20 个。

#function definition:

def Lorentzian(x, amp, cen, wid, n):
    f = 0
    for i in range( int(n) ):
        "lorentzian function: wid = half-width at half-max"
        f += (amp[i]/(1 + ((x-cen[i])/wid[i])**2))
    return f

#library import and model definition:

import lmfit as lf

lmodel = lf.Model(Lorentzian)

#The initial parameters for the model:
peaks_in_interval = np.array([2378, 2493, 2525, 2630, 2769])

number_of_peaks = len(peaks_in_interval)         
amplitude = width = np.zeros( number_of_peaks ) + 1
center = x[peaks_in_interval]

params = lmodel.make_params(x = x, amp = amplitude, cen = center, wid = width, n = number_of_peaks)

#This is the line that doesn't work:                        
result = lmodel.fit( y, params, x = x )

我已经开始尝试创建一个返回多洛伦兹，但我正在努力如何使其发挥作用......

我还发送 x、y 数组的数据。

DataSet for x and y

This is what the DataSet of x and y looks like.

Answer 1

您应该能够使用内置模型并使用前缀，如 manual 中所述。 。此外，最近在 mailinglist 上有一个关于非常相似主题的讨论。 .

您可以执行如下所示的操作。它还不太适合最后一个峰值，但您可能可以稍微调整一下起始值等。此外，由于您的基线并不完全平坦，当您使用 LinearModel 而不是 ConstantModel 时，它可能会有所改善，但我还没有尝试过。

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()