多指数衰减的 lmfit 模型

Question

我正在尝试使用 lmfit 包对我的数据进行拟合。但是我找不到任何用于多指数衰减的内置模型。我试图创建自己的函数，然后适应它。 我的代码如下：

import os
import time
import numpy as np
import pandas as pd
import lmfit
from lmfit.models import ExponentialModel, LinearModel
from lmfit import Model, Parameter, report_fit

def MultiExpDecay(tiempo,C1,tau1,C2,tau2):
    return C1*np.exp(-tiempo/tau1)+C2*np.exp(-tiempo/tau2)

def MultiExpDecay_fit():
    C1s = []
    C1s_error = []
    C2s = []
    C2s_error = []
    tau1s = []
    tau1s_error = []
    tau2s = []
    tau2s_error = []
    Fit_MultiExpDecays = []
    model = Model(MultiExpDecay, independent_vars=['tiempo'])
    for c in range(len(V_APD.columns)):
        xdat = tiempo.iloc[:, c]
        ydat = V_APD.iloc[:, c]
        pars = model.guess(ydat, x=xdat)
        fit = model.fit(ydat, pars, x=xdat)
        fit_values = model.eval(pars, x=xdat)
        Fit_MultiExpDecays.append(fit_values)
        for key in fit.params:
            if key == 'C1':
                #print(key, "=", out.params[key].value, "+/-", out.params[key].stderr)
                C1s.append(fit.params[key].value)
                C1s_error.append(fit.params[key].stderr)
            elif key == 'C2':
                C2s.append(fit.params[key].value)
                C2s_error.append(fit.params[key].stderr)
            elif key == 'tau1':
                tau1s.append(fit.params[key].value)
                tau1s_error.append(fit.params[key].stderr)
            elif key == 'tau2':
                tau2s.append(fit.params[key].value)
                tau2s_error.append(fit.params[key].stderr)
    Fit_MultiExpDecays = np.transpose(pd.DataFrame(Fit_MultiExpDecays, index=labels))
    C1 = np.transpose(pd.DataFrame(C1s, index = labels, columns = ['C1']))
    C2 = np.transpose(pd.DataFrame(C2s, index=labels, columns=['C2']))
    tau1 = np.transpose(pd.DataFrame(tau1s, index=labels, columns=['tau1']))
    tau2 = np.transpose(pd.DataFrame(tau2s, index=labels, columns=['tau2']))
    C1_error = np.transpose(pd.DataFrame(C1s_error, index = labels, columns=['C1 error']))
    C2_error = np.transpose(pd.DataFrame(C2s_error, index=labels, columns=['C2 error']))
    tau1_error = np.transpose(pd.DataFrame(tau1s_error, index=labels, columns=['tau1 error']))
    tau2_error = np.transpose(pd.DataFrame(tau2s_error, index=labels, columns=['tau2 error']))
    C1 = pd.concat([C1, C1_error])
    C2 = pd.concat([C2, C2_error])
    tau1 = pd.concat([tau1, tau1_error])
    tau2 = pd.concat([tau2, tau2_error])
    return C1, C2, tau1, tau2, Fit_MultiExpDecays

C1, C2, tau1, tau2, Fit_MultiExpDecays = MultiExpDEcay_fit():

出现错误但无法识别问题。

File "C:\ProgramData\Anaconda3\Lib\site-packages\lmfit\model.py", line 737, in guess
    raise NotImplementedError(msg)
NotImplementedError: guess() not implemented for Model

Answer 1

我支持@mikuszefski 的回答。你也可以使用两个lmfit.ExponentialModel，就像

import numpy as np
from lmfit.models import ExponentialModel

def dblexp( x, c1, l1, c2, l2 ):
    return c1 * np.exp( -x / l1 ) + c2 * np.exp( -x / l2 )
    
xl = np.linspace(0, 10, 101)
yl = dblexp(xl, 32.44, 0.81, 11.51, 9.22 ) + np.random.normal(0, 0.2, xl.size)

mymodel = ExponentialModel(prefix='e1_') + ExponentialModel(prefix='e2_') 

params = mymodel.make_params(e1_amplitude=10, e1_decay=10,
                             e2_amplitude=25, e2_decay=0.50)

result = mymodel.fit(yl, params, x=xl)
print( result.fit_report() )

但@mikuszefski 的回答或多或少是一样的。

对于实际问题：正如异常所说，Model 不知道如何猜测任意用户定义模型函数的合理参数值 - Model.guess() 方法默认引发此异常，必须由每个子类。 lmfit 为所提供的子模型（包括 ExponentialModel）提供了这种方法，但这些都是在理解这些特定模型的情况下编码的。

最后：拟合双指数衰减是出了名的容易出错，并且对于最小二乘方法（或至少 Levenberg-Marquardt 实现）来说是困难的。对于更实际的情况，我建议将数据的对数与双指数衰减的对数拟合。

Answer 2

是我错了还是你让你的生活太复杂了。虽然我自己没有使用lmfit，但我认为它的设计使您不必在此处完成所有编程。我认为它应该看起来像：

import matplotlib.pyplot as plt
from numpy import linspace
from numpy import fromiter
from numpy import exp
from numpy.random import normal
from lmfit import Model


def dblexp( x, c1, l1, c2, l2 ):
    return c1 * exp( -x / l1 ) + c2 * exp( -x / l2 )
    
xl = linspace(0, 10, 101)
yl = dblexp( xl, 32.44, 0.81, 11.51, 9.22 ) + normal(0, 0.2, xl.size)


mymodel = Model( dblexp ) 

### need some good guesses to start with
params = mymodel.make_params(c1=30, l1=1, c2=5, l2 =10)

result = mymodel.fit(yl, params, x=xl)
print( result.fit_report() )

fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.scatter( xl, yl )
ax.plot( xl, result.best_fit, 'r-' )
plt.show()

提供：

[[Model]]
    Model(dblexp)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 26
    # data points      = 101
    # variables        = 4
    chi-square         = 4.77426620
    reduced chi-square = 0.04921924
    Akaike info crit   = -300.239903
    Bayesian info crit = -289.779421
[[Variables]]
    c1:  32.5481660 +/- 0.18540661 (0.57%) (init = 30)
    l1:  0.79183180 +/- 0.00931390 (1.18%) (init = 1)
    c2:  11.6241194 +/- 0.16669056 (1.43%) (init = 5)
    l2:  9.11790608 +/- 0.19623957 (2.15%) (init = 10)
[[Correlations]] (unreported correlations are < 0.100)
    C(c2, l2) = -0.949
    C(l1, c2) = -0.822
    C(l1, l2) =  0.733
    C(c1, c2) = -0.652
    C(c1, l2) =  0.645
    C(c1, l1) =  0.251

和