# 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. – joriki Sep 9 '15 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? – wdkrnls Dec 15 '17 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$. – wdkrnls Dec 15 '17 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. – JJacquelin Dec 16 '17 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 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$. – JJacquelin Jan 21 at 8:03

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 ! – Claude Leibovici Jan 20 at 14:11
• @Claude Leibovici. Hi Claude ! Thank you for the compliment and for your interest in this subject. – JJacquelin Jan 20 at 16:57
• Thank you very much. Will try it shortly. – DrM Jan 21 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 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 at 14:32

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] = 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] = 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[1])

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] = 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] = 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] = 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[1])

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! – Ken Grimes Feb 15 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 at 19:55

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.