Python 计算随机分布的概率

Python 计算随机分布的概率，有了不完全伽马函数gamma和完全伽马函数Gamma_Half，就可以计算 x^2 的CDF值了。该值表示一个给定值是随机还是具有某种相关性。

函数本身十分短小：

def cdf(x: Union[Fraction, float], k: int) -> Fraction:
    """X2 cumulative distribution function.
    :param x: X2 value, sum (obs[i]-exp[i])**2/exp[i]
        for parallel sequences of observed and expected values.
    :param k: degrees of freedom >= 1; often len(data)-1
    """
    return (
        1 -
        gamma(Fraction(k, 2), Fraction(x/2)) /
        Gamma_Half(Fraction(k, 2))
    )

该函数包含了一些docstring注释用于解释参数。我们通过自由度和卡方值 x 创建了合适的Fraction对象。参数 x 可以是float值或者Fraction对象，这种灵活性可用于匹配那些完全使用浮点近似的示例。

可以使用Fraction(x/2).limit_denominator(1000)将x/2Fraction方法的大小限制为相对少量的数字。这会计算出一个正确的CDF值，而不是包含几十位数的庞大分数。

以下是从 x^2 表中调取的一些示例数据。更多相关信息，可参考维基百科词条chi-squared distribution。

执行以下命令以计算正确的CDF值：

>>> round(float(cdf(0.004, 1)), 2)
0.95
>>> cdf(0.004, 1).limit_denominator(100)
Fraction(94, 99)
>>> round(float(cdf(10.83, 1)), 3)
0.001
>>> cdf(10.83, 1).limit_denominator(1000)
Fraction(1, 1000)
>>> round(float(cdf(3.94, 10)), 2)
0.95
>>> cdf(3.94, 10).limit_denominator(100)
Fraction(19, 20)
>>> round(float(cdf(29.59, 10)), 3)
0.001
>>> cdf(29.59, 10).limit_denominator(10000)
Fraction(8, 8005)

在给定 x^2 和多个自由度的情况下，CDF函数生成的值与广泛使用的表格中的值一致。第一个示例展示了1个自由度时 x^2 为0.004的概率。第二个示例显示了1个自由度时 x^2 为10.38的概率。小的 x^2 值意味着预期结果和观测结果几乎没有差别。

下面是 x^2 表格中的一整行，由一个简单的生成器表达式计算得到：

>>> chi2 = [0.004, 0.02, 0.06, 0.15, 0.46, 1.07, 1.64, ...2.71, 3.84, 6.64, 10.83]
>>> act = [round(float(x), 3)
...    for x in map(cdf, chi2, [1]*len(chi2))]
>>> act
[0.95, 0.888, 0.806, 0.699, 0.498, 0.301, 0.2, 0.1, 0.05, 0.01, 0.001]

这些值显示了在1个自由度下，给定 x^2 和结果之间的相对似然度。与已发布的结果相比，计算结果在第三位小数上有一些细微差异。这意味着可以使用CDF计算结果代替在标准统计参考中查找到的 x^2 值。

函数CDF()给出了随机得到的 x^2 的概率值。

从已发布的表中可以看出，6个自由度下概率0.05的 x^2 值为12.5916。该CDF()函数的输出如下，显示出与已发布结果较为一致：

>>> round(float(cdf(12.5916, 6)), 2)
0.05

回顾先前示例，当时算出 x^2 的实际值是19.18。该值为随机的概率是：

>>> round(float(cdf(19.18, 6)), 5)
0.00387

当分母限制为1000时，概率为3/775，说明这些数据不太可能是随机的。这意味着我们可以拒绝零假设，并通过更多的分析来确定造成差异的可能原因。