using undetermined number of parameters in scipy function curve_fit

Question

First question: I'm trying to fit experimental datas with function of the following form:

f(x) = m_o*(1-exp(-t_o*x)) + ... + m_j*(1-exp(-t_j*x))

Currently, I don't find a way to have an undetermined number of parameters m_j, t_j, I'm forced to do something like this:

def fitting_function(x, m_1, t_1, m_2, t_2):
    return m_1*(1.-numpy.exp(-t_1*x)) + m_2*(1.-numpy.exp(-t_2*x)) 

parameters, covariance = curve_fit(fitting_function, xExp, yExp, maxfev = 100000)

(xExp and yExp are my experimental points)

Is there a way to write my fitting function like this:

def fitting_function(x, li):
    res = 0.
    for el in range(len(li) / 2):
        res += li[2*idx]*(1-numpy.exp(-li[2*idx+1]*x))
    return res

where li is the list of fitting parameters and then do a curve_fitting? I don't know how to tell to curve_fitting what is the number of fitting parameters. When I try this kind of form for fitting_function, I have errors like "ValueError: Unable to determine number of fit parameters."

Second question: Is there any way to force my fitting parameters to be positive?

Any help appreciated :)

Answer 1

See my question and answer here. I've also made a minimal working example demonstrating how it could be done for your application. I make no claims that this is the best way - I am muddling through all this myself, so any critiques or simplifications are appreciated.

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as pl

def wrapper(x, *args): #take a list of arguments and break it down into two lists for the fit function to understand
    N = len(args)/2
    amplitudes = list(args[0:N])
    timeconstants = list(args[N:2*N])
    return fit_func(x, amplitudes, timeconstants)


def fit_func(x, amplitudes, timeconstants): #the actual fit function
    fit = np.zeros(len(x))
    for m,t in zip(amplitudes, timeconstants):
        fit += m*(1.0-np.exp(-t*x))
    return fit

def gen_data(x, amplitudes, timeconstants, noise=0.1): #generate some fake data
    y = np.zeros(len(x))
    for m,t in zip(amplitudes, timeconstants):
        y += m*(1.0-np.exp(-t*x))
    if noise:
        y += np.random.normal(0, noise, size=len(x))
    return y


def main():
    x = np.arange(0,100)
    amplitudes = [1, 2, 3]
    timeconstants = [0.5, 0.2, 0.1]
    y = gen_data(x, amplitudes, timeconstants, noise=0.01)

    p0 = [1, 2, 3, 0.5, 0.2, 0.1]
    popt, pcov = curve_fit(lambda x, *p0: wrapper(x, *p0), x, y, p0=p0) #call with lambda function
    yfit = gen_data(x, popt[0:3], popt[3:6], noise=0)
    pl.plot(x,y,x,yfit)
    pl.show()
    print popt
    print pcov

if __name__=="__main__":
    main()

A word of warning, though. A linear sum of exponentials is going to make the fit EXTREMELY sensitive to any noise, particularly for a large number of parameters. You can test that by adding even a small amount of noise to the data generated in the script - even small deviations cause it to get the wrong answer entirely while the fit still looks perfectly valid by eye (test with noise=0, 0.01, and 0.1). Be very careful interpreting your results even if the fit looks good. It's also a form that allows for variable swapping: the best fit solution is the same even if you swap any pairs of (m_i, t_i) with (m_j, t_j), meaning your chi-square has multiple identical local minima that might mean your variables get swapped around during fitting, depending on your initial conditions. This is unlikely to be a numeriaclly robust way to extract these parameters.

To your second question, yes, you can, by defining your exponentials like so:

m_0**2*(1.0-np.exp(-t_0**2*x)+...

Basically, square them all in your actual fit function, fit them, and then square the results (which could be negative or positive) to get your actual parameters. You can also define variables to be between a certain range by using different proxy forms.