I have written a function to compute the Laplace transform of a function using scipy.integrate.quad. It is not a very sophisticated function and currently performs poorly on the probability density function of an Erlang distribution.
I have included all my work below. I first compute the Laplace transform and then the inverse in order to compare it to the original p.d.f. of the Erlang. I use mpmath for this. The mpmath.invertlaplace is not the problem as it manages to convert the closed-form Laplace transform back to the original p.d.f. quite perfectly.
Please help me understand what the problem is with my numerical laplace transform. I get the following error but have not been able to resolve it.
IntegrationWarning: The occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be underestimated. a=0,b=np.inf,limit=limit)[0]
PLOT
CODE
import numpy as np
import matplotlib.pyplot as plt
import math as m
import mpmath as mp
import scipy.stats as st
from scipy.integrate import quad
def get_laplace(func,limit=10000):
'''
Returns laplace transfrom function
'''
def laplace(s):
'''Numerical laplace transform'''
# Seperate into real and imaginary parts
x = np.real(s)
y = np.imag(s)
def real_func(t):
return m.exp(-x*t)*m.cos(y*t)*func(t)
def imag_func(t):
return m.exp(-x*t)*m.sin(y*t)*func(t)
real_integral = quad(real_func,
a=0,b=np.inf,limit=limit)[0]
imag_intergal = quad(real_func,
a=0,b=np.inf,limit=limit)[0]
return complex(real_integral,-imag_intergal)
return laplace
def L_erlang(s,lam,k):
'''
Closed form laplace transform of Erlang or Gamma distribution.
'''
return (lam/(lam+s))**k
if __name__ == '__main__':
# Setup Erlang probability density function
k = 5
lam = 1
pdf = st.erlang(a=k,scale=1/lam).pdf
# Laplace transforms
Lnum = get_laplace(pdf) # numerical approximation
L = lambda s: L_erlang(s,lam,k) # closed form
# Use mpmath library to perform inverse laplace
# Invserse transfrom on numerical laplace function
invLnum = lambda t: mp.invertlaplace(Lnum,t,
method='dehoog',
dps=5,
degree=3)
# Invserse transfrom on closed-form laplace function
invL = lambda t: mp.invertlaplace(L,t,
method='dehoog',
dps=5,
degree=3)
# Grid to visualise
T = np.linspace(0.1,10,10)
# Get values of inverse transforms
lnum = np.array([invLnum(t) for t in T])
l = np.array([invL(t) for t in T])
# Plot
plt.plot(T,lnum,label='from numerical laplace')
plt.plot(T,l,label='from closed-form laplace')
plt.plot(T,pdf(T),label='original pdf',linestyle='--')
plt.legend(loc='best')
plt.show()
UPDATE
After two cups of VERY strong coffee, I managed to see the obvious mistake and make the code work. It quite embarrassing actually. Have a look at this line of code:
imag_intergal = quad(real_func,
a=0,b=np.inf,limit=limit)[0]
Hmmm, real_func hey? So it should read:
imag_intergal = quad(imag_func_func,
a=0,b=np.inf,limit=limit)[0]
The one gets this lovely plot:
Conclusion
So why go through all the trouble to numerically perform the Laplace transform of something that we have a closed form solution for. That is because the interest lies somewhere else. We do not have a closed expression for the conditional future lifetime distribution which is similar to the hazard function. Lets call it h. Then for the Erlang distribution erl = st.erlang(a=k,scale=1/lam) that has been active for tau units of time we have h = lambda t: erl.pdf(t+tau)/erl.sf(tau). This distribution can be used as a holding time in a Semi-Markov Model (SMP). To analyse the transient behaviour of the SMP one use Laplace transforms. Usually only pdfs are used but now I can use hazard functions. Its pretty cool because it means one can model transient behaviour without assuming everything is new.
