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