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.
