python - scipy stats binom cdf 返回 nan

Question

如果我理解正确，scipy.stats 离散分布的 cdf 应该返回给定参数值的概率总和。

因此，scipy.stats.binom(7000000000, 0.5).cdf(6999999999) 应该返回几乎恰好为 1 的值，因为在 70 亿次试验中，有 50/50 的概率在 70 亿减 1 或更少的情况下取得成功是非常确定的。相反，我得到了 np.nan。事实上，对于提供给 .cdf 的任何值，除了 70 亿本身(或更多)，我都会返回 np.nan。

这是怎么回事？ scipy.stats 发行版可以处理的数字是否有一些限制，但文档中没有？

Answer 1

TL;博士

内部计算时缺乏浮点精度。虽然 scipy 是一个 Python 库，但它的核心是用 C 语言编写的，并使用 C 数字类型。

---

举个例子:

import scipy.stats

for i in range (13):
    trials = 10 ** i
    print(f"i: {i}\tprobability: {scipy.stats.binom(trials, 0.5).cdf(trials - 1)}")

输出是:

i: 0    probability: 0.5
i: 1    probability: 0.9990234375
i: 2    probability: 0.9999999999999999
i: 3    probability: 0.9999999999999999
i: 4    probability: 0.9999999999999999
i: 5    probability: 0.9999999999999999
i: 6    probability: 0.9999999999999999
i: 7    probability: 0.9999999999999999
i: 8    probability: 0.9999999999999999
i: 9    probability: 0.9999999999999999
i: 10   probability: nan
i: 11   probability: nan
i: 12   probability: nan

原因在于二项分布的 CDF 公式(我无法嵌入图像，所以这里是维基链接:https://en.wikipedia.org/wiki/Binomial_distribution

在 scipy 源代码中，我们会看到对此实现的引用:http://www.netlib.org/cephes/doubldoc.html#bdtr

在它的深处涉及到按trials划分(incbet.c, line 375: ai = 1.0/a; 这里叫做a,但是 nwm)。如果您的 trials 太大，这个除法的结果会很小，以至于当我们将这个小数字加到另一个不是那么小的数字时，它实际上并没有改变，因为我们缺乏浮点精度这里(目前只有 64 位)。然后，经过更多的算术运算，我们尝试从一个数字中得到对数，但它等于零，因为它没有在应该改变的时候改变。而log(0)没有定义，等于np.nan。