Python 通过高斯和多指数衰减的卷积进行曲线拟合的代码

Question

我正在开发一个代码，用于将数据与一个模型进行拟合，该模型是两个函数的卷积（具有多指数衰减的高斯函数 exp(Ax)+exp(Bx)+...）。基本上只有高斯 and/or 高斯修改 https://en.wikipedia.org/wiki/Exponentially_modified_Gaussian_distribution 的拟合在 Lmfit 中工作得很好但是使用内置卷积（即如果使用两个函数的 np.convolve Lmfit 不起作用。

我在互联网上尝试了很多例子，到目前为止我意识到我的函数 returns inf 或 nan 值以及数据在卷积中的使用不是等间距的。我通过使用卷积的数学表达式和使用 scipy.optimize.curve_fit 找到了解决这个问题的弯路。但这是一个非常笨拙和耗时的方法，我想找到一种方法通过使用使它更复杂和通用两个函数的卷积并使用 lmfit，我可以更轻松地控制参数。

该数据集也包含在评论中，供您参考。

w=0.1 # is constant
def CONVSum(x,w,*p):
    n=np.int(len(p)/3)
    A=p[:n]
    B=p[n:2*n]
    C=p[2*n:3*n]
# =======================================================================
#     below formula is derived as mathematical expression of convoluted multi exponential components with a gaussian distribution based on the instruction given in http://www.np.ph.bham.ac.uk/research_resources/programs/halflife/gauss_exp_conv.pdf
# ======================================================================
    fnct=sum(np.float64([A[i]*np.exp(-B[i]*((x-C[i])-(0.5*np.square(w)*B[i])))*(1+scipy.special.erf(((x-C[i])-(np.square(w)*B[i]))/(np.sqrt(2)*w))) for i in range(n)]))
    fnct[np.isnan(fnct)]=0
    fnct[fnct<1e-12]=0
    return fnct
N=4 #number of exponential functions to be fitted
params = np.linspace(1, 0.0001, N*3); #parameters for a multiple exponential                                                                                                                        
popt,pcov = curve_fit(CONVSum,x,y,p0=params,
    bounds=((0,0,0,0,-np.inf,-np.inf,-np.inf,-np.inf,-3,-3,-3,-3),
           (1,1,1,1, np.inf, np.inf, np.inf, np.inf, 3, 3, 3, 3)),
           maxfev = 1000000)

fitted data with curve fitt

非常感谢任何有关高斯卷积和多重指数衰减拟合的帮助或提示，我更喜欢使用 lmfit，因为我可以很好地识别参数并将它们相互关联。

理想情况下，我想让我的数据与参数相匹配，其中一些参数在数据集之间共享，一些被延迟 (+off_set)。

Answer 1

好吧，你的脚本有点难以阅读，并且有很多与你的问题无关的内容。您的 exgauss 函数无法防止无穷大。 np.exp(x) for x>~ 710 会给出Inf，拟合将无法进行。

Answer 2

这里是有问题的固化适配代码的等价物。我设法通过使用 and here 中非常棒的说明和信息来创建它。但是还是有待开发。

    # =============================================================================
    #     below formula is drived as mathematical expresion of convoluted multi exponential components with a gausian distribution based on the instruction given in http://www.np.ph.bham.ac.uk/research_resources/programs/halflife/gauss_exp_conv.pdf
    # =============================================================================

def CONVSum(x,params):
    fnct=sum(
            np.float64([
                    (params['amp%s_%s'%(n,i)].value)*np.exp(-(params['dec%s_%s'%(n,i)].value)*((x-(params['cen%s_%s'%(n,i)].value))-
                     (0.5*np.square((params['sig%s_%s'%(n,i)].value))*(params['dec%s_%s'%(n,i)].value))))*
                    (1+scipy.special.erf(((x-(params['cen%s_%s'%(n,i)].value))-(np.square((params['sig%s_%s'%(n,i)].value))*
                    (params['dec%s_%s'%(n,i)].value)))/(np.sqrt(2)*(params['sig%s_%s'%(n,i)].value)))) for n in range(N) for i in wav 
                    ])
            )
    fnct=fnct/fnct.max()
    return fnct
    # =============================================================================
    #  this global fit were adapted from  

    # it is of very important thet we can identify the shared parameteres for datasets
    # =============================================================================


def objective(params, x, data):
    """ calculate total residual for fits to several data sets"""
    ndata = data.shape[0]
    resid = 0.0*data[:]
    # make residual per data set
    resid = data- CONVSum(x,params)
    # now flatten this to a 1D array, as minimize() needs
    return resid.flatten()


# selec datasets
x  = df[949].index
data =df[949].values


# create required sets of parameters, one per data set
N=4 #number of exponential decays
wav=[949]  #the desired data to be fitted
fit_params = Parameters()
for i in wav:
        for n in range(N):
            fit_params.add( 'amp%s_%s'%(n,i), value=1, min=0.0,  max=1)
            fit_params.add( 'dec%s_%s'%(n,i), value=0.5, min=-1e10,  max=1e10)
            fit_params.add( 'cen%s_%s'%(n,i), value=0.1, min=-3.0,  max=1000)
            fit_params.add( 'sig%s_%s'%(n,i), value=0.1, min=0.05, max=0.5)

# now we constrain some values to have the same value
# for example assigning sig_2, sig_3, .. sig_5 to be equal to sig_1
for i in wav:
        for n in (1,2,3):
            print(n,i)
            fit_params['sig%s_%s'%(n,i)].expr='sig0_949'
            fit_params['cen%s_%s'%(n,i)].expr='cen0_949'

# it will run the global fit to all the data sets
result = minimize(objective, fit_params, args=(x,data)) 
report_fit(result.params)

# plot the data sets and fits
plt.close('all')
plt.figure()
for i in wav:
    y_fit = CONVSum(x,result.params)
    plt.plot(x, data, 'o-', x, y_fit, '-')
    plt.xscale('symlog') 
plt.show()

fitted data with convolution of multi exponential and gausian

不幸的是，拟合结果不是很令人满意，我仍在寻找一些改进的建议。

Python 通过高斯和多指数衰减的卷积进行曲线拟合的代码

Python code for curve fitting by convolution of a gausian and multi exponential decay

python

gaussian

curve-fitting

convolution

exponential