Python中使用数值积分的拉普拉斯变换精度非常差

Question

我编写了一个函数来计算函数的拉普拉斯变换scipy.integrate.quad。它不是一个非常复杂的函数，目前在 Erlang 分布的概率密度函数上表现不佳。

我在下面列出了我的所有工作。我首先计算拉普拉斯变换，然后计算逆变换，以便将其与 Erlang 的原始 pdf 进行比较。我用mpmath这个。这mpmath.invertlaplace不是问题，因为它成功地将封闭式拉普拉斯变换非常完美地转换回原始 pdf。

请帮助我理解我的数值拉普拉斯变换的问题是什么。我收到以下错误，但无法解决它。

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]

嗯，real_func嘿？所以它应该是：

imag_intergal = quad(imag_func_func,
                     a=0,b=np.inf,limit=limit)[0]

人们得到了这个可爱的情节：

结论 那么，为什么要费尽心思对我们已经有闭合形式解的事物进行拉普拉斯变换呢？那是因为兴趣在其他地方。我们没有类似于危险函数的条件未来寿命分布的封闭表达式。就这样吧h。erl = st.erlang(a=k,scale=1/lam)那么对于已经活跃了tau一定时间单位的Erlang 发行版，我们有h = lambda t: erl.pdf(t+tau)/erl.sf(tau)。该分布可用作半马尔可夫模型 (SMP) 中的保持时间。为了分析 SMP 的瞬态行为，可以使用拉普拉斯变换。通常只使用 pdf，但现在我可以使用危险函数。它非常酷，因为这意味着人们可以对瞬态行为进行建模，而无需假设一切都是新的。

Answer 1

再会。我重复原始问题的更新部分，因为这是问题的解决方案。这样，问题就可以标记为已解决。

更新

喝了两杯非常浓的咖啡后，我设法看到了明显的错误并使代码正常工作。其实也挺尴尬的。看看这行代码：

imag_intergal = quad(real_func,
                     a=0,b=np.inf,limit=limit)[0]

嗯，real_func嘿？所以它应该是：

imag_intergal = quad(imag_func_func,
                     a=0,b=np.inf,limit=limit)[0]

人们得到了这个可爱的情节：

结论 那么，为什么要费尽心思对我们已经有闭合形式解的事物进行拉普拉斯变换呢？那是因为兴趣在其他地方。我们没有类似于危险函数的条件未来寿命分布的封闭表达式。就这样吧h。erl = st.erlang(a=k,scale=1/lam)那么对于已经活跃了tau一定时间单位的Erlang 发行版，我们有h = lambda t: erl.pdf(t+tau)/erl.sf(tau)。该分布可用作半马尔可夫模型 (SMP) 中的保持时间。为了分析 SMP 的瞬态行为，可以使用拉普拉斯变换。通常只使用 pdf，但现在我可以使用危险函数。它非常酷，因为这意味着人们可以对瞬态行为进行建模，而无需假设一切都是新的。