# Fit sum of exponentials

I am looking for an algorithm that fits a sum of exponential. for example I have something like this: $$y(x)=ae^{−bx}+c+de^{-fx}+h$$ and I want to find a,b,c,d,f and h values. Of this sum I have only the x and y values that belong to the curve represented by the previous model. Is there any paper that explains that?

• For a start, you should combine $c$ and $h$ into a single additive constant; you can't determine them separately. Sep 9, 2015 at 20:40

A direct method of fitting (no initial guess, no iteration) for the function : $$y(x)=a+be^{px}+ce^{qx}$$ is summarized below (parameters to be computed : $a,b,c,p,q$ ). It works as well in case of negative $p, q$ : Instead of minimizing the absolute deviations, the variant below minimizes the relatives deviations : The theory of this method is given in the paper :https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales (in French).

The method for the function $$y(x)=a+be^{px}+ce^{qx}+de^{rx}$$ is also available, but not published yet. Contact the author if interrested.

• Hello @JJacquelin. I would be interested in seeing your tri-exponential regression result. I tried to follow your paper, to implement my own, but I'm getting stuck with all imaginary solutions to p, q, r. Have you by any chanced published it in the interim? Dec 15, 2017 at 23:07
• The integral equation I came up with was: $y = A(SSS)+B(-SS)+C(S)+Dx^2+Ex+F$ where $S$* refers to the number of integrations. This lead me to solve the equations $A=p q r, B=qr+pr+pq, C=p+q+r$. Dec 15, 2017 at 23:08
• @ wdkrnls : Your calculus is correct. Well done ! So, the linear regression leads to approximates of $A,B,C,D,E,F$. Solving the cubic polynomial equation for $p,q,r$ should give the three real roots. But if the scatter on the data is too large and/or if the number of points is not large enough, the numerical calculus fails. Number of points and scatter are very critical factors because the successive numerical integrations are likely to introduce too much deviations. Dec 16, 2017 at 7:36
• @JJacquelin Hi. Did you publish an explicit solution for the 3 exponent case? And have you posted any implementations? Thank you. Will cite your work when we publish, of course.
– DrM
Jan 19, 2020 at 15:31
• @wdkrnls I suppose that your calculus is for the case where $a=0$. If $a\neq 0$ there should be a $x^3$ term in the integral equation. A cause of imaginary solutions to $p,q,r$ might be the inversion of $A=p+q+r$ and $C=pqr$. Jan 21, 2020 at 8:03

Edit: some people have asked me about the license of the code below. The license is MIT as specfified in this Github repository.

I would like to add another solution, after the useful discussion with @JJacquelin. I reworked my previous solution based on linear algebra and derivatives, to work with integrals instead. It also works for $$n$$ exponentials.

Assuming the data set $$(y, x)$$;

Compute the $$n$$ integrals of $$y$$ ; $$\int y, \int^2{y}, \cdots , , \int^{n-1}{y} , \int^n{y}$$.

Compute the $$n-1$$ powers of the independent variable $$x$$ ; $$x, x^2, \cdots , x^{n-1}$$.

Solve linear least squares problem for $$[a_{1}, a_{2}, \cdots, a_{n-1}, a_{n}, k_{n-1}, k_{n-2} , \cdots , k_{2}, k_{1}, k_{0}]$$ in:

$$y = \left[ \int y , \int^2{y} , \cdots, \int^{n-1}{y}, \int^n{y} , x^{n-1}, x^{n-2}, \cdots , x^{2}, x, 1 \right] \left[ \begin{array}{c} a_{1} \\ a_{2} \\ \vdots \\ a_{n-1} \\ a_{n} \\ k_{n-1} \\ k_{n-2} \\ \vdots \\ k_{2} \\ k_{1} \\ k_{0} \end{array} \right]$$

$$y = Y A \\ Y^T Y A = Y^T y \\ A = (Y^T Y)^{-1}Y^T y$$

Then use the first $$n$$ elements of $$A$$ to form the transfer matrix of the integral equation:

$$\bar{A} = \left[ \begin{array}{cccccc} a_1 & a_2 & a_3 & \cdots & a_{n-1} & a_{n} \\ 1 & 0 & 0 & \cdots & 0 & 0 \\ 0 & 1 & 0 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 1 & 0 \end{array} \right]$$

Compute the eigenvalues of $$\bar{A}$$ to obtain the solution of the characteristic polynomial.

The eigenvalues $$\lambda _n$$ are then the values of the exponentials in:

$$y = p_1 e^{\lambda _1 x} + p_2 e^{\lambda _2 x} + \cdots + p_{n-1} e^{\lambda _{n-1} x} + p_n e^{\lambda _n x}$$

To obtain the missing $$p_n$$, calculate numerically the exponential terms $$e^{\lambda _n}$$ and solve the least squares problem for $$[p_1, p_2, \cdots, p_{n-1}, p_n ]$$ in

$$y = [e^{\lambda _1 x}, e^{\lambda _2 x} , \cdots, e^{\lambda _{n-1} x}, e^{\lambda _n x}] \left[ \begin{array}{c} p_1 \\ p_2 \\ \vdots \\ p_{n-1} \\ p_n \end{array} \right]$$

$$y = X P \\ X^T X P = X^T y \\ P = (X^T X)^{-1}X^T y$$

Now you have your $$n$$ exponential model parameters $$p_n$$ and $$\lambda _n$$.

This solution based on integrals, looks much more similar to @JJacquelin's solutions for cases $$n = 2$$ and $$n = 3$$, but now generalized for any $$n$$.

Here is the Matlab code for a test case with $$n = 4$$ (you can test it online here):

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 = eig([A(1), A(2), A(3), A(4); 1, 0, 0, 0; 0, 1, 0, 0; 0, 0, 1, 0]);
lambdas
%lambdas =
%  -2.9991
%  -1.9997
%   0.5000
%   0.1500

% 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


## Variations

If it is desired to fit a sum of exponentials plus constant $$k$$:

$$y = k + p_1 e^{\lambda _1 x} + p_2 e^{\lambda _2 x} + \cdots + p_{n-1} e^{\lambda _{n-1} x} + p_n e^{\lambda _n x}$$

Just add $$x^n$$ to the first least squares problem and a constant to the second least squares problem:

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

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

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

lambdas = eig([A(1), A(2), A(3); 1, 0, 0; 0, 1, 0]);
lambdas
%lambdas =
%  -2.9991
%  -1.9997
%   0.5000

% get exponentials multipliers
X = [ones(size(x)), exp(lambdas(1)*x), exp(lambdas(2)*x), exp(lambdas(3)*x)];
P = pinv(X)*y;
P
%P =
%  -0.9996
%   4.0043
%   1.9955
%   4.9999


### Sine and Cosine

If it is desired to fit a sum of exponentials plus a $$cosine$$ or a $$sine$$, the algorithm is exactly the same. Sines and cosines manifest themselves as a pair of complex conjugate eigen-values of $$\hat{A}$$.

For example, if we fit one exponential and one cosine:

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

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

% get exponentials lambdas
Y = [iy1, iy2, iy3, x.^2, x, ones(size(x))];
A = pinv(Y)*y;
%lambdas =
%   0.5000 +      0i
%  -0.0000 + 0.1500i
%  -0.0000 - 0.1500i

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

% get exponentials multipliers
X = [exp(lambdas(1)*x), exp(lambdas(2)*x), exp(lambdas(3)*x)];
P = pinv(X)*y;
P
%P =
%   5.0001 + 0.0000i
%   2.0000 + 0.0001i
%   2.0000 - 0.0001i


Note the cosine multiplier $$4$$ appears in $$P$$ as $$2$$ times $$2$$.

If we were to fit a $$sine$$ instead we would get:

y  =  5*exp(0.5*x) + 4*sin(0.15*x);

% other stuff

lambdas =

0.5000 +      0i
-0.0000 + 0.1500i
-0.0000 - 0.1500i

P =

5.0001e+00 + 5.2915e-16i
-4.6543e-05 - 1.9999e+00i
-4.6543e-05 + 1.9999e+00i


Note now the $$2$$ times $$2$$ appears in the complex part of $$P$$.

### Exponential multiplied by $$x$$

If it is desired to fit a sum of exponentials, some with repeated $$\lambda$$ but multiplied by $$x$$, again the exact same method can be used. This scenario manifest itself as repeated eigen-values of $$\hat{A}$$:

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

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

% get exponentials lambdas
Y = [iy1, iy2, iy3, x.^2, x, ones(size(x))];
A = pinv(Y)*y;

lambdas = eig([A(1), A(2), A(3); 1, 0, 0; 0, 1, 0]);
lambdas
%lambdas =
%  -2.9991
%  -2.9991 NOTE : repeated means x * e^{lambda * x}
%   0.5000

% get exponentials multipliers
X = [exp(lambdas(1)*x), x.*exp(lambdas(2)*x), exp(lambdas(3)*x)];
P = pinv(X)*y;
P
%P =
%   4.0001
%   1.9955
%   5.0000


## What if I don't know how many exponentials fit my data?

If the structure of the model to fit is unknown, one can simply run the algorithm for a large $$n$$ and compute the singular values of $$Y$$. By looking at the singular values, we can infer how many of them are actually needed to fit most of the data. From this we can select a minimum $$n$$ and based in the singular values of $$\hat{A}$$ we can also deduce the structure of the model (how many exponentials, sines, cosines, if there are exponentials multiplied times $$x$$, etc.).

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

% calculate n integrals of y and n-1 powers of x
n = 10;
iy = zeros(length(x), n);
xp = zeros(length(x), n-1);
iy(:,1) = cumtrapz(x, y);
xp(:,1) = ones(size(x));
for ii=2:1:n
iy(:, ii) = cumtrapz(x, iy(:, ii-1));
xp(:, ii) = xp(:, ii-1) .* x;
end

% calculate singular values of Y
Y = [iy, xp];
ysv = svd(Y);

% scale singular values to percentage
ysv_scaled = 100 * ysv ./ sum(ysv)
bar(ysv_scaled);
%ysv_scaled =
%   6.9606e+01
%   2.3696e+01
%   4.2812e+00
%   1.5760e+00
%   6.8682e-01
%   1.3106e-01
%   2.0743e-02
%   2.5832e-03
%   2.5537e-04
%   ... even smaller ones

% select n principal components above a threshold of 0.1%
thres = 0.1;
n = sum(ysv_scaled > 0.1) / 2
% n = 3 (6 singular values) covers about 99.97% of all components contributions
covers = sum(ysv_scaled(1:2*n))
% covers = 99.976

• Bravo ! Smart theoretical work. Nevertheless, in practical applications I am afraid that successive numerical integrations will give more and more inconsistant results. But this is much better than succesive numerical differentiations due to the noise. Aug 30, 2020 at 13:59
• Note : To be more complet the function to fit should be the sum of $n$ exponentials plus one unknown constant. But this doesn't change the method. Only a few changes in the linear algorithm. Best regards. Aug 30, 2020 at 14:18
• Yes, just added it. Actually the method works for sines, cosines, and exponentials multiplied by powers of $x$. It depends on the eigenvalues of $\hat{A}$, I just don't have the time to put up all variations. Actually by analyzing the principal components (singular value decomposition) of $Y^T Y$ with a large $n = m$ where $m >> n$ we could infer how many exponentials are needed to get a useful fit with minimum $n$. Aug 30, 2020 at 15:08
• Very well. The method works for many other functions: Weibull, logistic, gaussian, mixed polynomial and sinusoidal, etc. A lot of examples are shown in the paper fr.scribd.com/doc/14674814/Regressions-et-equations-integrales . Aug 30, 2020 at 15:27
• Mon francais ne marche pas tres bien, mais je vais essayer! Aug 30, 2020 at 15:29

In response to questions raised in comments about the three exponents case, the method of regression based on integral equation is roughly explain below, with a numerical example.  • Always impressive ! Jan 20, 2020 at 14:11
• @Claude Leibovici. Hi Claude ! Thank you for the compliment and for your interest in this subject. Jan 20, 2020 at 16:57
• Thank you very much. Will try it shortly.
– DrM
Jan 21, 2020 at 14:25
• Aside, I am looking into transforms that provide the rate constants as poles. There is the Pady-Laplace of course. The fourier transform of x e^{-r x} and variations, seem interesting also. Have you worked with methods of that sort? What are your thoughts? Thank you
– DrM
Jan 21, 2020 at 14:30
• Ojne of the physically most interesting questions in a large class of problems, is how many exponential rates are there? It is about how many states are there in the system, to which the number of rates is the answer.
– DrM
Jan 21, 2020 at 14:32

The solution is on wikipedia.

If we want to fit to n exponentials, follow the steps:

Calculate $$n$$ derivatives of $$y$$ ; $$y^{(n)}, y^{(n-1)}, \cdots , \ddot y , \dot y , y$$

Solve linear least squares problem for $$[-a_{n-1}, \cdots, -a_2, -a_1, -a_0 ]$$ in

$$y^{(n)} = [y^{(n-1)} , \cdots, \ddot y, \dot y, y] \left[ \begin{array}{c} -a_{n-1} \\ \vdots \\ -a_2 \\ -a_1 \\ -a_0 \end{array} \right]$$

$$y^{(n)} = Y A \\ Y^T Y A = Y^T y^{(n)} \\ A = (Y^T Y)^{-1}Y^T y^{(n)}$$

Obtain the roots of the characteristic polynomial (using Matlab roots or similar):

$$z^n + a_{n-1}z^{n-1} + \cdots + a_2 z^2 + a_1z + a_0$$

The roots $$\lambda _n$$ are then the values of the exponentials in

$$y = p_1 e^{\lambda _1 x} + p_2 e^{\lambda _2 x} + \cdots + p_{n-1} e^{\lambda _{n-1} x} + p_n e^{\lambda _n x}$$

To obtain the missing $$p_n$$, calculate numerically the exponential terms $$e^{\lambda _n}$$ and solve the least squares problem for $$[p_1, p_2, \cdots, p_{n-1}, p_n ]$$ in

$$y = [e^{\lambda _1 x}, e^{\lambda _2 x} , \cdots, e^{\lambda _{n-1} x}, e^{\lambda _n x}] \left[ \begin{array}{c} p_1 \\ p_2 \\ \vdots \\ p_{n-1} \\ p_n \end{array} \right]$$

$$y = X P \\ X^T X P = X^T y \\ P = (X^T X)^{-1}X^T y$$

Now you have your $$n$$ exponential model parameters $$p_n$$ and $$\lambda _n$$.

I believe this approach is similar of not the same as the one proposed by @JJacquelin, but you have a generic algorithm for $$n$$ exponentials.

If your $$y$$ data is noise, I would recommend applying a non-causal low pass filter first.

Here is the Matlab code (you can test it online here):

clear all;
clc;
% get data
dx = 0.01;
x  = (dx:dx:1.5)';
y  =  5*exp(0.5*x) + 4*exp(-3*x) + 2*exp(-2*x);
% calculate derivatives
dy1 = diff(y)/dx;
dy2 = diff(dy1)/dx;
dy3 = diff(dy2)/dx;
% fix derivatives lengths
dy1 = [dy1(1)-(dy1(2)-dy1(1)); dy1];
dy2 = [dy2(1)-2*(dy2(2)-dy2(1)); dy2(1)- (dy2(2)-dy2(1)); dy2];
dy3 = [dy3(1)-3*(dy3(2)-dy3(1)); dy3(1)-2*(dy3(2)-dy3(1)); dy3(1)- (dy3(2)-dy3(1)); dy3];
% get exponentials lambdas
Y = [dy2, dy1, y];
A = pinv(Y)*dy3;
lambdas = roots([1; -A]);
lambdas
%lambdas =
%  -3.0336
%  -1.9933
%   0.4989
% get exponentials multipliers
X = [exp(lambdas(1)*x), exp(lambdas(2)*x), exp(lambdas(3)*x)];
P = pinv(X)*y;
P
%P =
%   3.9461
%   2.0514
%   5.0067


Here is are implementations of the two and three exponent cases, from Jean Jacquelin's work.

#!/usr/bin/python3

'''
Implementation of two and three exponent cases from Jean Jacquelin's
REGRESSION et EQUATION INTEGRALE, available at
https://www.scribd.com/doc/14674814/Regressions-et-equations-integrales
'''

__author__    = "M Nelson"
__date__      = "January 19, 2020"
__version__   = "0.1"

import sys
import argparse

import numpy as np
from math import sqrt
from scipy.linalg import lstsq
from scipy.optimize import curve_fit

def biexpfit( xdata, ydata, debug=False ):

x = np.array(xdata)
y = np.array(ydata)
if debug:
print( 'x', x )
print( 'y', y )

S = np.empty_like(y)
S = 0
S[1:] = np.cumsum(0.5 * (y[1:] + y[:-1]) * np.diff(x))
if debug:
print('S', S )

SS = np.empty_like(y)
SS = 0
SS[1:] = np.cumsum(0.5 * (S[1:] + S[:-1]) * np.diff(x))
if debug:
print('SS', SS )

x2 = x * x
x3 = x2 * x
x4 = x2 * x2

M = [ [sum(SS*SS),  sum(SS*S), sum( SS*x2 ), sum(SS*x), sum(SS)],
[sum(SS*S),   sum(S*S),  sum(S*x2),    sum(S*x), sum(S) ],
[sum(SS*x2),  sum(S*x2), sum(x4),      sum(x3),  sum(x2) ],
[sum(SS*x),   sum(S*x),  sum(x3),      sum(x2),  sum(x) ],
[sum(SS),     sum(S),    sum(x2),      sum(x),   len(xdata) ] ]

if debug:
print( 'M' )
for n in range(5):
print( M[n] )

Ycol = np.array( [ sum(SS*y), sum(S*y), sum(x2*y), sum(x*y), sum(y) ] )

(A,B,C,D,E),residues,rank,singulars = list( lstsq( M, Ycol ) )
if debug:
print( 'A-E', A, B, C, D, E )

'''
Minv = np.linalg.inv(M)
A,B,C,D,E = list( np.matmul(Minv,Ycol) )
'''

p = (1/2.)*(B + sqrt(B*B+4*A))
q = (1/2.)*(B - sqrt(B*B+4*A))
if debug:
print( 'p,q', p, q )

beta = np.exp(p*x)
eta = np.exp(q*x)

betaeta = beta * eta

L = [ [ len(xdata), sum(beta), sum(eta) ],
[ sum(beta),  sum(beta*beta), sum(betaeta) ],
[ sum(eta),   sum(betaeta), sum(eta*eta)] ]

Ycol = np.array( [ sum(y), sum(beta*y), sum(eta*y) ] )

(a,b,c),residues,rank,singulars = list( lstsq( L, Ycol ) )
if debug:
print( 'a,b,c', a,b,c )

'''
Linv = np.linalg.inv(L)
a,b,c = list( np.matmul( Linv, Ycol ) )
'''

# sort in ascending order (fastest negative rate first)
(b,p),(c,q) = sorted( [[b,p],[c,q]], key=lambda x: x)

return a,b,c,p,q

def triexpfit( xdata, ydata, debug=False ):

x = np.array(xdata)
y = np.array(ydata)
if debug:
print( 'x', x )
print( 'y', y )

S = np.empty_like(y)
S = 0
S[1:] = np.cumsum(0.5 * (y[1:] + y[:-1]) * np.diff(x))
if debug:
print('S', S )

SS = np.empty_like(y)
SS = 0
SS[1:] = np.cumsum(0.5 * (S[1:] + S[:-1]) * np.diff(x))
if debug:
print('SS', SS )

SSS = np.empty_like(y)
SSS = 0
SSS[1:] = np.cumsum(0.5 * (SS[1:] + SS[:-1]) * np.diff(x))
if debug:
print('SSS', SSS )

x2 = x * x
x3 = x2 * x
x4 = x3 * x
x5 = x4 * x
x6 = x5 * x

M = [ [ sum(SSS*SSS), sum(SSS*SS), sum(SSS*S), sum(SSS*x3), sum(SSS*x2), sum(SSS*x), sum(SSS) ],
[ sum(SSS*SS), sum(SS*SS), sum(SS*S), sum(SS*x3), sum(SS*x2), sum(SS*x), sum(SS) ],
[ sum(SSS*S), sum(SS*S), sum(S*S), sum(S*x3), sum(S*x2), sum(S*x), sum(S) ],
[ sum(SSS*x3), sum(SS*x3), sum(S*x3), sum(x6), sum(x5), sum(x4), sum(x3) ],
[ sum(SSS*x2), sum(SS*x2), sum(S*x2), sum(x5), sum(x4), sum(x3), sum(x2) ],
[ sum(SSS*x), sum(SS*x), sum(S*x), sum(x4), sum(x3), sum(x2), sum(x) ],
[ sum(SSS), sum(SS), sum(S), sum(x3), sum(x2), sum(x), len(xdata) ] ]

if debug:
print( 'M' )
for n in range(7):
print( M[n] )

Ycol = [ sum(y*SSS), sum(y*SS), sum(y*S), sum(y*x3), sum(y*x2), sum(y*x), sum(y) ]

if debug:
print( 'Y', Ycol )

'''
Minv = np.linalg.inv(M)
A,B,C,D,E,F,G = list( np.matmul( Minv, Ycol ) )
print( 'A-G linalg', A, B, C, D, E, F, G )
'''

(A,B,C,D,E,F,G),residues,rank,singulars = list( lstsq( M, Ycol ) )
if debug:
print( 'A-G', A, B, C, D, E, F, G )

p,q,r = np.roots( [1.,-C,-B,-A] )
p,q,r = sorted( [p,q,r] )
if debug:
print( 'p,q,r', p, q, r )

L = [ [ len(xdata), sum(np.exp(p*x)),sum(np.exp(q*x)),sum(np.exp(r*x)) ],
[ sum(np.exp(p*x)), sum(np.exp(2*p*x)),sum(np.exp((p+q)*x)),sum(np.exp((p+r)*x)) ],
[ sum(np.exp(q*x)), sum(np.exp((p+q)*x)),sum(np.exp(2*q*x)),sum(np.exp((q+r)*x)) ],
[ sum(np.exp(r*x)), sum(np.exp((p+r)*x)),sum(np.exp((q+r)*x)),sum(np.exp(2*r*x)) ] ]

Ycol = [ sum(y), sum(y*np.exp(p*x)),sum(y*np.exp(q*x)),sum(y*np.exp(r*x)) ]

(a,b,c,d),residues,rank,singulars = list( lstsq( L, Ycol ) )
if debug:
print( 'a,b,c,d', a,b,c,d )

# sort in ascending order (fastest negative rate first)
(b,p),(c,q),(d,r) = sorted( [[b,p],[c,q],[d,r]], key=lambda x: x)

return a,b,c,d, p,q,r

def __regression_tests__(debug=False):

x = np.linspace( 0, 2., 100 )

a, b = 1., 2.
p = -5.
print( 'inputs', a, b, p )

y = a + b * np.exp( p*x )

a, b, p = expfit( x, y )
print( 'fits', a, b, p )

a, b, c = 1.0, .5, 0.25
p,q = -10., -5.
print( 'inputs', a,b,c,p,q )

y = a + b * np.exp( p*x ) + c * np.exp( q* x )

a,b,c,p,q =  biexpfit( x, y )
print( 'fits', a,b,c,p,q )

a, b, c, d = -1., 5., 4., 2.
p, q, r = .5, -3., -2.
print( 'inputs', a,b,c,d, p,q,r )

x = np.linspace( .1, 1.5, 15 )
print( x )

y = a + b * np.exp( p*x ) + c * np.exp( q* x ) + d * np.exp( r* x )
print( y )

a,b,c,d, p,q,r = triexpfit( x, y, debug )
print( 'fits', a,b,c,d, p,q,r )

• Wow, amazing work. Thanks a lot! Feb 15, 2020 at 12:46
• @KenGrimes Thank you, but really this is Jean's solution. I just coded it. Also helpful were comments from another contributor. A caveat to note is that while the 2 term case is pretty good, and especially useful for the baseline, the 3 term case is not always better than the usual minimization methods.
– DrM
Feb 17, 2020 at 19:55

In order to answer to the comment from Yanqi Huang the method of fitting $$y=be^{px}+ce^{qx}$$ is shown below. This is the same as in the first answer but with one line and column suppressed related to the suppression of the parameter $$a$$. There is no longer $$a\:x^2$$ in the linear integral equation. Fall on your discussion while looking at how to fit Weibull distribution (PDF or CDF). As you may know, these days we look at many exponentials curves which commonly result in the sum of local exponentials. I am no expert but I was giving it a try. I try solving a sum of Weibull distributions (which gives exponentials when the shape parameter equals one). I did not find immediately a minimizing solution, so I was just looking to your method in case. And because a few equations are worth many words: If any idea or suggestion, gladly appreciated :)

In all cases, thanks so much for this post. I guess I will spend the day implementing your approach.