使用 python 和 lmfit 测量 irf 的迭代再卷积拟合

Question

我正在尝试使用 python 和 lmfit.

将带卷积的指数衰减函数拟合到测得的仪器响应

我是 python 的新手，我正在尝试遵循 https://groups.google.com/group/lmfit-py/attach/90f51c25ebb39a52/deconvol_exp2.py?part=0.1&authuser=0 中的代码。

import numpy as np
from lmfit import Model 
import matplotlib.pyplot as plt
import requests

# Load data
url = requests.get('https://groups.google.com/group/lmfit-py/attach/73a983d40ad945b1/tcspcdatashifted.csv?part=0.1&authuser=0')

x,decay1,irf=np.loadtxt(url.iter_lines(),delimiter=',',unpack=True,dtype='float')

plt.figure(1)
plt.semilogy(x,decay1,x,irf)
plt.show()

# Define weights
wWeights=1/np.sqrt(decay1+1)


# define the double exponential model
def jumpexpmodel(x,tau1,ampl1,tau2,ampl2,y0,x0,args=(irf)):
    ymodel=np.zeros(x.size) 
    t=x
    c=x0
    n=len(irf)
    irf_s1=np.remainder(np.remainder(t-np.floor(c)-1, n)+n,n)
    irf_s11=(1-c+np.floor(c))*irf[irf_s1.astype(int)]
    irf_s2=np.remainder(np.remainder(t-np.ceil(c)-1,n)+n,n)
    irf_s22=(c-np.floor(c))*irf[irf_s2.astype(int)]
    irf_shift=irf_s11+irf_s22   
    irf_reshaped_norm=irf_shift/sum(irf_shift)    
    ymodel = ampl1*np.exp(-(x)/tau1)
    ymodel+= ampl2*np.exp(-(x)/tau2)  
    z=Convol(ymodel,irf_reshaped_norm)
    z+=y0
    return z

# convolution using fft (x and h of equal length)
def Convol(x,h):
    X=np.fft.fft(x)
    H=np.fft.fft(h)
    xch=np.real(np.fft.ifft(X*H))
    return xch    

# assign the model for fitting    
mod=Model(jumpexpmodel)

在定义拟合的初始参数时，我收到错误。

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

#initialize the parameters - showing error
pars = mod.make_params(tau1=10,ampl1=1000,tau2=10,ampl2=1000,y0=0,x0=10,args=irf)
pars['x0'].vary =True
pars['y0'].vary =True

print(pars)

# fit this model with weights, initial parameters
result = mod.fit(decay1,params=pars,weights=wWeights,method='leastsq',x=x)

# print results
print(result.fit_report())

# plot results
plt.figure(5)
plt.subplot(2,1,1)
plt.semilogy(x,decay1,'r-',x,result.best_fit,'b')
plt.subplot(2,1,2)
plt.plot(x,result.residual)
plt.show()

根据 lmfit.model 的文档，我怀疑这是因为参数 irf 在模型中如何定义为 args=(irf)。我试图将 irf 传递给 model 而不是 params。我也尝试过使用 **kwargs.

将irf合并到model中进行卷积并拟合数据的正确方法是什么？

Answer 1

我相信您想将 irf 视为模型函数的附加 自变量 - 您传递给函数但未处理的值作为拟合中的变量。

为此，只需将模型函数的签名 jumpexpmodel() 修改为更简单的

def jumpexpmodel(x, tau1, ampl1, tau2, ampl2, y0, x0, irf):

函数体很好（实际上 args=(irf) 不会起作用，因为您需要解包 args -- 这里的签名确实是您想要的）。

然后告诉lmfit.Model()，irf是一个自变量——默认是第一个参数是唯一的自变量：

mod = Model(jumpexpmodel, independent_vars=('x', 'irf'))

那么在做参数的时候，不要包含irf或者args:

pars = mod.make_params(tau1=10, ampl1=1000, tau2=10, ampl2=1000, y0=0, x0=10)

而是现在将 irf 与 x 一起传递给 mod.fit():

result = mod.fit(decay1, params=pars, weights=wWeights, method='leastsq', x=x, irf=irf)

您的程序的其余部分看起来不错，结果拟合效果相当好，给出了

的报告

[[Model]]
    Model(jumpexpmodel)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 138
    # data points      = 2797
    # variables        = 6
    chi-square         = 3795.52585
    reduced chi-square = 1.35991610
    Akaike info crit   = 865.855713
    Bayesian info crit = 901.473529
[[Variables]]
    tau1:   50.4330421 +/- 0.68246203 (1.35%) (init = 10)
    ampl1:  2630.30664 +/- 20.1552948 (0.77%) (init = 1000)
    tau2:   225.392872 +/- 2.75674753 (1.22%) (init = 10)
    ampl2:  523.257894 +/- 12.4451921 (2.38%) (init = 1000)
    y0:     20.7975212 +/- 0.14165429 (0.68%) (init = 0)
    x0:    -9.70588133 +/- 0.12597936 (1.30%) (init = 10)
[[Correlations]] (unreported correlations are < 0.100)
    C(tau2, ampl2) = -0.947
    C(tau1, ampl2) = -0.805
    C(tau1, tau2)  =  0.706
    C(tau1, x0)    = -0.562
    C(ampl1, x0)   =  0.514
    C(tau1, ampl1) = -0.453
    C(tau2, y0)    = -0.426
    C(ampl2, y0)   =  0.314
    C(ampl2, x0)   =  0.291
    C(tau2, x0)    = -0.260
    C(tau1, y0)    = -0.212
    C(ampl1, tau2) =  0.119

和这样的情节：

使用 python 和 lmfit 测量 irf 的迭代再卷积拟合

Iterative reconvolution fitting with measured irf using python and lmfit

python

convolution

non-linear-regression

lmfit

deconvolution