lmfit model for multi exponential decay

Question

I am trying to perform a fitting for my data using the lmfit package. However I couldnt find any built-in model for a multi-exponential decay. I have tried to create my own function, and then to fit it. My code is the following:

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

An error is raised but cannot identify the problem.

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 is designed such that you do not have to do all the programming you have done here. it should just look like:

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

Providing:

[[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

and

Answer 2

You could also use two lmfit.ExponentialModel, as with

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

For the actual question: As the exception says, Model does not know how to guess reasonable parameter values for an arbitrary, user-defined model function - the Model.guess() method raises this exception by default and must be overwritten by each subclass. lmfit provides this method for the provided submodels (including ExponentialModel), but these are each coded with an understanding of those particular models.

Finally: fitting double-exponential decay is notoriously error-prone and difficult for least-squares approaches (or at least Levenberg-Marquardt implementations). For more realistic cases, I would recommend fitting the log of the data with the log of the double-exponential decay.