Неправильное кумулятивное распределение, полученное при интегрировании PDF - PullRequest
Новая
       92

Неправильное кумулятивное распределение, полученное при интегрировании PDF

0 голосов
/ 06 апреля 2020

У меня есть PDF (код ниже), и я пытаюсь получить его кумулятивная функция распределения (CDF). Я могу использовать стандартный подход интеграции PDF с 0. (это минимально допустимое значение) для увеличения x значений. Это работает, как и ожидалось:

enter image description here

, где PDF отображается красным, а "пошаговый" CDF синим. Но это оставляет меня с таблицей (x, y) данных, и я хочу фактическую функцию , которая описывает этот CDF. Поэтому я перехожу к интеграции PDF от 0. до переменной y, чтобы получить это выражение. Я делаю это с WolframAlpha (WA) и нахожу:

enter image description here

Я пишу эту функцию в свой код, используя верхняя неполная гамма-функция , и я получаю это:

enter image description here

, где WA интегрированный CDF выделен оранжевым цветом. Я попробовал нижнюю неполную гамма-функцию, но результаты хуже.

Я почти уверен, что WA CDF написан без ошибок, поэтому я не уверен, что я делаю здесь неправильно. 

import numpy as np
from scipy.special import gamma, gammaincc
from scipy import integrate
import matplotlib.pyplot as plt

def main():
    xx = np.linspace(0., 5., 100)
    yy = PDF(xx)
    plt.plot(xx, yy, c='r')

    # Find stepwise CDF
    cummul = []
    for x in xx:
        cummul.append([x, integrate.quad(PDF, 0., x)[0]])
    plt.plot(*np.array(cummul).T)

    # WolframAlpha's integral
    y = CDF(xx)
    plt.plot(xx, y)


def PDF(y):
    a = (343. / 15.) * np.sqrt(7. / (2. * np.pi))
    b, c = 5. / 2., -7. / 2.
    return a * (y ** b) * np.exp(c * y)


def CDF(x):
    """
    From WolframAlpha
    """
    a = (343. / 15.) * np.sqrt(7. / (2. * np.pi))
    b, c = 5. / 2., -7. / 2.
    return a * (
        x**b * (-c * x)**(-b) * (gammaincc(b + 1, -c * x) - b * gamma(b))) / c

if __name__ == '__main__':
    main()

1 Ответ

1 голос
/ 06 апреля 2020

scipy.special.gammaincc - регуляризованная верхняя неполная гамма-функция . Если неполная гамма-функция - это гамма (a, x), то переквалифицированная неполная гамма-функция - это гамма (a, x) / gamma (a). Неполная гамма-функция в формуле WA не является регуляризованной версией. Вы получите ожидаемый результат, если умножите gammaincc(b + 1, -c * x) на gamma(b + 1). То есть, измените оператор возврата на 

    return a * (
        x**b * (-c * x)**(-b) * (gamma(b + 1) * gammaincc(b + 1, -c * x) - b * gamma(b))) / c

, который также может быть записан 

    return a * (
        x**b * (-c * x)**(-b) * gamma(b + 1) * (gammaincc(b + 1, -c * x) - 1)) / c

Если я использую последнюю версию, и изменим график интеграла Wolfram Alpha на 

    plt.plot(xx, y, '--', linewidth=4, alpha=0.6)

Я получаю этот участок:

plot

...