多指数衰减的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

您也可以使用两个

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

对于实际问题：正如异常所说，

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

和