我编写了一个函数,用scipy.integrate.quad
计算函数的拉普拉斯变换。它不是一个非常复杂的函数,目前在Erlang分布的概率密度函数上表现不佳
我已经把我所有的工作都包括在下面了。我首先计算拉普拉斯变换,然后求逆,以便将其与Erlang的原始p.d.f.进行比较。我使用mpmath
来实现这一点。mpmath.invertlaplace
不是问题所在,因为它能够非常完美地将闭合形式的拉普拉斯变换转换回原始的p.d.f
请帮助我理解数值拉普拉斯变换的问题所在。我收到以下错误,但无法解决它
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]
绘图
代码
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()
更新
喝了两杯浓咖啡后,我终于发现了明显的错误,并使代码正常工作。事实上这很尴尬。请看一下这行代码:
imag_intergal = quad(real_func,
a=0,b=np.inf,limit=limit)[0]
嗯,嘿?因此,它应该是:
imag_intergal = quad(imag_func_func,
a=0,b=np.inf,limit=limit)[0]
其中一个得到了这个可爱的情节:
结论
那么,为什么要费尽心机,对我们有封闭形式解的东西进行拉普拉斯变换呢。这是因为利益在其他地方。我们没有一个类似于危险函数的条件未来寿命分布的闭合表达式。我们把它叫做h
。然后对于活动了tau
个时间单位的Erlang分布erl = st.erlang(a=k,scale=1/lam)
,我们有h = lambda t: erl.pdf(t+tau)/erl.sf(tau)
。这种分布可以用作半马尔可夫模型(SMP)中的保持时间。要分析SMP的瞬态行为,请使用拉普拉斯变换。通常只使用PDF,但现在我可以使用危险函数。这很酷,因为这意味着人们可以在不假设一切都是新的情况下模拟瞬时行为
目前没有回答
相关问题 更多 >
编程相关推荐