binom cdf返回nan

Question

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

因此，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；DR

在内部计算过程中缺乏浮点精度。虽然C是一个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：distribution的链接)。

在this源中，我们可以看到对这个实现的修改：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。