python - 通过集成 PDF 获得的错误累积分布

Question

我有一个 PDF（下面的代码），我正在尝试获取它的累积分布函数(CDF)。我可以使用标准方法将 PDF 从 0.（它是允许的最小值）集成到递增x值。这按预期工作：

其中 PDF 为红色，“逐步”CDF 为蓝色。但这给我留下了一张数据表，我想要 描述这个 CDF(x, y)的实际函数。所以我转向将 PDF 从 0. 整合到一个变量y来获得这个表达式。我用 WolframAlpha (WA) 做这个，我发现：

我使用上面的不完全伽马函数将此函数写入我的代码中，我得到了这个：

其中 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()

标签： pythonscipycdf

Answer 1

scipy.special.gammaincc是正则化的上不完全伽马函数。如果不完全 gamma 函数是 gamma(a, x)，则正则化的不完全 gamma 函数是 gamma(a, x)/gamma(a)。WA 公式中的不完全 gamma 函数不是正则化版本。如果乘以 ，您将得到预期的gammaincc(b + 1, -c * x)结果gamma(b + 1)。也就是把return语句改成

    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)

我得到这个情节：