多指数衰减的 lmfit 模型

lmfit model for multi exponential decay

我正在尝试使用 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

是我错了还是你让你的生活太复杂了。虽然我自己不使用 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

我支持@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 实现)来说很困难。对于更实际的情况,我建议将数据的日志与双指数衰减的日志进行拟合。