Exponential fit - sum of exponentials

Question

I want to fit the following data:

x(time) y(data)    
0.75;   19.33    
1;      19.04   
1.25;   17.21   
1.5;    12.98   
1.75;   11.59   
2;  9.26   
2.25;   7.66   
2.5;    6.59    
2.75;   5.68    
3;  5.1    
3.25;   4.36   
3.5;    4.43     
3.75;   3.58    
4;  3.01    
4.25;   3.24     
4.5;    3.58     
4.75;   3.13     
5;  3.88     
5.25;   3.19     
5.5;    3.58     
5.75;   3.64   

using the following code:

#read text file
data = pd.read_table('episode_5_prova.txt', sep='\t')
#DataFrame
df = pd.DataFrame(data)
#Define your function
def func(x, a, b, c, d, e):
return a*np.exp(-b*x) + c*np.exp(-d*x) + e

#convert dataframe into numpy array
df0=df['time']
x=df0.as_matrix()
df1=df['bi']
y=df1.as_matrix()

# Using the python scipy function curve_fit with function and input variables
popt, pcov = curve_fit(func, x, y)
a, b, c, d, e= popt
fit = func(x, a, b, c, d, e)

fig, ax = plt.subplots()
ax.plot(x, fit, color='r', lw=3)
ax.plot(x, y,'g.')
observed_values=scipy.array(y)
expected_values=scipy.array(fit)
plt.xlim(0,25)
plt.ylim(0,20)

print(a,b,c,d,e)

print(scipy.stats.chisquare(observed_values, f_exp=expected_values, ddof=3))
plt.show()

I obtain the following plot: first fit However, for the purpose of my work, I need to fix my parameters b and c as: b=0.000431062, d=0.000580525 but i don't obtained a good fit as follows: second fit

Does anyone have a suggestion? Thanks

Answer 1

This is a duplicate of this question.

In short, you need an eigenvalue solver, a least squares solver and follow the steps outlined in the solution.

There is Matlab code available that you can use as an example to implement in other languages (e.g. Python).

clear all;
clc;
% get data
dx = 0.02;
x  = (dx:dx:1.5)';
y  =  5*exp(0.5*x) + 4*exp(-3*x) + 2*exp(-2*x) - 3*exp(0.15*x);

% calculate integrals
iy1 = cumtrapz(x, y);
iy2 = cumtrapz(x, iy1);
iy3 = cumtrapz(x, iy2);
iy4 = cumtrapz(x, iy3);

% get exponentials lambdas
Y = [iy1, iy2, iy3, iy4, x.^3, x.^2, x, ones(size(x))];
A = pinv(Y)*y;
%lambdas =
%  -2.9991
%  -1.9997
%   0.5000
%   0.1500

lambdas = eig([A(1), A(2), A(3), A(4); 1, 0, 0, 0; 0, 1, 0, 0; 0, 0, 1, 0]);
lambdas

% get exponentials multipliers
X = [exp(lambdas(1)*x), exp(lambdas(2)*x), exp(lambdas(3)*x), exp(lambdas(4)*x)];
P = pinv(X)*y;
P
%P =
%   4.0042
%   1.9955
%   4.9998
%  -2.9996