Обобщенное распределение вероятностей Student-T, которое я кодировал в Python, не интегрируется с 1 (в некоторых случаях)

Question

Я пытался реализовать асимметричное обобщенное распределение t в Python, чтобы смоделировать некоторую финансовую отдачу. Я основал свой код на формулах, найденных в Википедии , и использовал бета-дистрибутив от scipy.

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. Во втором случае, однако, я не понимаю, почему значение так высоко.

Кто-нибудь знает, в чем проблема?

Это графики нормального (синего) и плохого (оранжевого)

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,999999999999999978 и 0,99752026308335162. Причина, по которой они не равны 1, заключается в том, что CDF определяется как интеграл от -infinity до бесконечности PDF.