【问题标题】:Use Python lmfit with a variable number of parameters in function在函数中使用带有可变数量参数的 Python lmfit
【发布时间】:2019-03-02 20:15:42
【问题描述】:

我正在尝试将复杂的气相色谱信号反卷积为单独的高斯信号。这是一个示例,其中虚线表示我尝试去卷积的信号。

我能够使用 scipy.optimize.curve_fit 编写代码来执行此操作;然而,一旦应用于真实数据,结果就不可靠了。我相信能够为我的参数设置界限会改善我的结果,所以我正在尝试使用 lmfit,它允许这样做。我在让 lmfit 使用可变数量的参数时遇到问题。我正在使用的信号可能具有任意数量的基础高斯分量,因此我需要的参数数量会有所不同。我在这里找到了一些提示,但仍然无法弄清楚......

Creating a python lmfit Model with arbitrary number of parameters

这是我目前正在使用的代码。代码将运行,但当模型拟合时,参数估计值不会改变。有谁知道我怎样才能让我的模型工作?

import numpy as np
from collections import OrderedDict
from scipy.stats import norm
from lmfit import Parameters, Model

def add_peaks(x_range, *pars):
    y = np.zeros(len(x_range))
    for i in np.arange(0, len(pars), 3):
        curve = norm.pdf(x_range, pars[i], pars[i+1]) * pars[i+2]
        y = y + curve
    return(y)

# generate some fake data
x_range = np.linspace(0, 100, 1000)
peaks = [50., 40., 60.]
a = norm.pdf(x_range, peaks[0], 5) * 2
b = norm.pdf(x_range, peaks[1], 1) * 0.1
c = norm.pdf(x_range, peaks[2], 1) * 0.1
fake = a + b + c

param_dict = OrderedDict()

for i in range(0, len(peaks)):
    param_dict['pk' + str(i)] = peaks[i]
    param_dict['wid' + str(i)] = 1.
    param_dict['mult' + str(i)] = 1.

# In case, you'd like to see the plot of fake data
#y = add_peaks(x_range, *param_dict.values())
#plt.plot(x_range, y)
#plt.show()

# Initialize the model and fit
pmodel = Model(add_peaks)

params = pmodel.make_params()
for i in param_dict.keys():
    params.add(i, value=param_dict[i])

result = pmodel.fit(fake, params=params, x_range=x_range)
print(result.fit_report())

【问题讨论】:

    标签: python signal-processing modeling lmfit deconvolution


    【解决方案1】:

    我认为您最好使用lmfits 构建复合模型的能力。
    也就是说,用

    定义的单个峰
    from scipy.stats import norm
    def peak(x, amp, center, sigma):
        return amp * norm.pdf(x, center, sigma)
    

    (另见lmfit.models.GaussianModel),您可以建立一个具有许多峰的模型:

    npeaks = 3
    model = Model(peak, prefix='p1_')
    for i in range(1, npeaks):
         model = model + Model(peak, prefix='p%d_' % (i+1))
    
    params = model.make_params()
    

    现在model 将是 3 个高斯函数的总和,为该模型创建的 params 将具有类似 p1_ampp1_centerp2_amp、...的名称,您可以添加 sensible初始值和/或界限和/或约束。

    根据您的示例数据,您可以将初始值传递给make_params,例如

    params = model.make_params(p1_amp=2.0, p1_center=50., p1_sigma=2, 
                               p2_amp=0.2, p2_center=40., p2_sigma=2, 
                               p3_amp=0.2, p3_center=60., p3_sigma=2)
    
    result = model.fit(fake, params, x=x_range)
    

    【讨论】:

      【解决方案2】:

      我在这里找到了解决方案:

      https://lmfit.github.io/lmfit-py/builtin_models.html#example-3-fitting-multiple-peaks-and-using-prefixes

      在上面的代码的基础上,以下完成了我想要做的事情......

      from lmfit.models import GaussianModel
      
      gauss1 = GaussianModel(prefix='g1_')
      gauss2 = GaussianModel(prefix='g2_')
      gauss3 = GaussianModel(prefix='g3_')
      gauss4 = GaussianModel(prefix='g4_')
      gauss5 = GaussianModel(prefix='g5_')
      
      gauss = [gauss1, gauss2, gauss3, gauss4, gauss5]
      prefixes = ['g1_', 'g2_', 'g3_', 'g4_', 'g5_']
      
      mod = np.sum(gauss[0:len(peaks)])
      pars = mod.make_params()
      
      for i, prefix in zip(range(0, len(peaks)), prefixes[0:len(peaks)]):
          pars[prefix + 'center'].set(peaks[i])
      
      init = mod.eval(pars, x=x_range)
      out = mod.fit(fake, pars, x=x_range)
      print(out.fit_report(min_correl=0.5))
      out.plot_fit()
      plt.show()
      

      【讨论】:

        猜你喜欢
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 2018-08-14
        • 1970-01-01
        • 1970-01-01
        • 2011-12-24
        • 2012-04-18
        相关资源
        最近更新 更多