python - 我用 Python 编写的广义 Student-T 概率分布没有积分为 1(在某些情况下)

Question

我一直在尝试在 Python 中实现偏斜广义 t 分布来模拟一些财务返回。我的代码基于在 Wikipedia 上找到的公式，我使用了 scipy 的 Beta 分布。

from scipy.special import beta
import numpy as np
from math import sqrt

def sgt(x, params):
# This function accepts an array of 5 parameters [mu, sigma, lambda, p, q]  
    mu, sigma, lam, p, q = params

    v = (q**(-1/p)) / (sqrt((3*lam*lam + 1)*beta(3/p, q-2/p)/beta(1/p, q) - 4*lam*lam*(beta(2/p, q-1/p)/(beta(1/p, q)))**2))
    m = 2*v*sigma*lam*q**(1/p)*beta(2/p, q - 1/p) / beta(1/p, q)
    fx = p / (2*v*sigma*(q**(1/p))*beta(1/p, q)*((abs(x-mu+m)**p/(q*(v*sigma)**p*(lam*np.sign(x-mu+m)+1)**p + 1)+1)**(1/p + q)))

    return fx

现在，该函数似乎对某些参数集工作得很好，但对其他参数集却很糟糕。

例如:

dx = 0.001
x_axis = np.arange(-10, 10, dx)

ok_parameters = [0, 2, 0, 3, 8]
bad_parameters = [0, 2, 0, 1.05, 2.1]

ok_distribution = sgt(x_axis, ok_parameters)
bad_distribution = sgt(x_axis, bad_parameters)

如果我尝试计算这两个数字的积分:

a = np.sum(ok_distribution*dx)
b = np.sum(bad_distribution*dx)

我得到结果 a = 1.0013233154393804 和 b = 2.2799746093533346。 现在，理论上这两个都应该为 1，但我假设因为我对积分进行了近似，所以该值不会总是正好为 1。但是在第二种情况下，我不明白为什么该值如此之高。

有人知道问题出在哪里吗？

These are the graphs of the ok distribution (blue) and bad distribution (orange)

Answer 1

我相信您的定义 sgt 中只是有一个拼写错误(虽然我找不到确切的位置)。这是一个有效的实现。

%matplotlib inline
import matplotlib.pyplot as plt
from scipy.special import beta
import numpy as np
from math import sqrt
from typing import Union
from scipy import integrate

# Generalised Student T probability Distribution
def generalized_student_t(x:Union[float, np.ndarray], mu:float, sigma:float, 
                          lam:float, p:float, q:float) \
        -> Union[float, np.ndarray]:


    v = q**(-1/p) * ((3*lam**2 + 1)*(beta(3/p, q - 2/p)/beta(1/p,q)) - 4*lam**2*(beta(2/p, q - 1/p)/beta(1/p,q))**2)**(-1/2)

    m = 2*v*sigma*lam*q**(1/p)*beta(2/p,q - 1/p)/beta(1/p,q)   

    fx = p  / (2*v*sigma*q**(1/p)*beta(1/p,q)*(abs(x-mu+m)**p/(q*(v*sigma)**p)*(lam*np.sign(x-mu+m)+1)**p + 1)**(1/p + q))

    return fx

def plot_cdf_pdf(x_axis:np.ndarray, pmf:np.ndarray) -> None:
    """
    Plot the PDF and CDF of the array returned from the function.
    """
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
    ax1.plot(x_axis, pmf)
    ax1.set_title('PDF')
    ax2.plot(x_axis, integrate.cumtrapz(x=x_axis, y=pmf, initial = 0))
    ax2.set_title('CDF')
    pass


dx = 0.0001
x_axis = np.arange(-10, 10, dx)

# Create the Two
distribution1 = generalized_student_t(x=x_axis, mu=0, sigma=1, lam=0, p=2, q=100)
distribution2 = generalized_student_t(x=x_axis, mu=0, sigma=2, lam=0, p=1.05, q=2.1)

plot_cdf_pdf(x_axis=x_axis, pmf=distribution1)
plot_cdf_pdf(x_axis=x_axis, pmf=distribution2)

我们还可以检查 PDF 的积分是否为 1

integrate.simps(x=x_axis, y = distribution1)
integrate.simps(x=x_axis, y = distribution2)

我们可以看到积分的结果是0.99999999999999978和0.99752026308335162。它们不完全为 1 的原因是 CDF 被定义为从 PDF 的 -infinity 到 infinity 的积分。