scipy stats binom cdf 返回 nan答案

【问题标题】：scipy stats binom cdf returns nanscipy stats binom cdf 返回 nan
【发布时间】：2018-11-08 08:02:19
【问题描述】：

如果我理解正确，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 发行版可以处理的数量是否存在文档中没有的限制？

【问题讨论】：

有趣。底层函数是scipy.special.bdtr。对于n=2**31-1，bdtr(n-1, n, 0.5) 返回 0.9999999999999999。对于n=2**31，它返回nan。
我怀疑这是 scipy 中的一个错误。
我在 github 上创建了一个问题：github.com/scipy/scipy/issues/9454。结果可能是重复的； scipy.stats.binom、scipy.special.bdtr 和 scipy.special.betainc 有几个相关问题。

标签： python scipy

【解决方案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 公式（我无法嵌入图像，所以这里是 wiki 的链接：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。

【讨论】：

我认为答案更简单：试验次数由bdtr 转换为 32 位整数。但是 10**10 溢出了一个 32 位整数，incbet 为负数 a 或 b 返回 NaN。尽管 incbet 确实存在精度问题 (github.com/scipy/scipy/issues/5503)，但它们在这里没有发挥作用。