替代允许 64 位或更多位的 numpy.random.binomial？答案

【问题标题】：Alternative to numpy.random.binomial that allows 64 bits or more?替代允许 64 位或更多位的 numpy.random.binomial？
【发布时间】：2014-05-20 17:10:17
【问题描述】：

我需要根据二项分布随机生成一系列数字。 Numpy 的随机套件提供了一种方法来执行此操作，但不幸的是，它似乎仅限于处理 n 值的 32 位数字，我想使用该范围之外的值。 64 位就足够了，尽管任意更高的精度也可以。

示例输出：

>>> np.random.binomial(1<<40, 0.5)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "mtrand.pyx", line 3506, in mtrand.RandomState.binomial (numpy\random\mtrand\mtrand.c:16555)
OverflowError: Python int too large to convert to C long

我可以使用其他替代方法吗？或者让 numpy 在这个随机生成器内部使用 64 位数字的方法？

或者我需要自己骑上马鞍并滚动吗？

（正如 Robert Klein 所指出的，numpy 确实在除 Windows 之外的 64 位平台上运行 64 位；不幸的是，我使用的是 Windows）。

【问题讨论】：

你能再解释一下吗？ np.random.binomial( np.iinfo(np.int64).max , 0.5) 在我的笔记本电脑上运行良好。

标签： python python-3.x numpy random

【解决方案1】：

在 C long 整数为 64 位的机器上，numpy.random.binomial() 将接受并生成 64 位整数。除 Windows 外，大多数 64 位平台都是如此。例如，在我的 64 位 OS X 机器上：

[~]
|11> np.random.binomial(1 << 40, 0.5)
549755265539

[~]
|12> np.random.binomial(1 << 40, 0.5) > (1<<32)
True

或者，如果您卡在 Windows 上，请考虑使用 Normal Approximation 进行二项分布。在这么大的n 上，近似值应该很好。

def approx_binomial(n, p, size=None):
    gaussian = np.random.normal(n*p, sqrt(n*p*(1-p)), size=size)
    # Add the continuity correction to sample at the midpoint of each integral bin.
    gaussian += 0.5
    if size is not None:
        binomial = gaussian.astype(np.int64)
    else:
        # scalar
        binomial = int(gaussian)
    return binomial

【讨论】：

有趣。不幸的是，我使用的是 Windows，所以它没有多大帮助。
考虑将Normal Approximation 用于二项分布。无论如何，它在这么大的n 上可能有更好的行为。 np.random.normal(n*p, np.sqrt(n*p*(1-p))).astype(np.int64)
您能否将其发布为答案，以便我接受？
将其添加到此答案中。谢谢！
我添加了一个连续性校正，发现我需要以正常方式转换为 int 而不是使用.astype，但否则效果很好。谢谢，罗伯特。

【解决方案2】：

有一个精确且有效的采样器来生成二项式 (n, 1/2) 随机变量，即使 n 很大，如 Bringmann 等人所述。 2014. 与该论文中描述的算法等效的算法如下：

如果 n 小于 4，则生成 n 个无偏随机位（零或一）并返回它们的总和。否则，如果 n 是奇数，则将 ret 设置为该算法的结果 n = n - 1，然后将一个无偏随机位的值添加到 ret，然后返回 ret。
将 m 设置为 floor(sqrt(n)) + 1。
（首先，从二项式曲线的包络中采样。）生成无偏随机位（零或一），直到以这种方式生成零。将 k 设置为以这种方式生成的数量。
将s设为[0,m)中的一个整数，均匀随机选择，然后将i设为k em>*m + s.
将 ret 设置为 n/2+i 或 n/2 - i - 1 等概率。
（二，接受或拒绝ret。）如果ret ret > n，去到第 3 步。
有概率选择(n, ret)*m*2^{(k - n)+2}，返回ret。否则，进入第3步。（这里，choose(n,k)是一个二项式系数。下面的Python代码计算了这个概率的对数。）

（请注意，Bringmann 论文的算法比这更复杂，部分是为了避免由于精度有限而导致的舍入误差，并且 Farach-Colton 和 Tsai [2015] 展示了如何解决抽样二项式的问题（n, p) 任何p 的随机变量都可以简化为对二项式(n, 1/2) 变量进行采样的问题。有关详细信息，请参阅那些论文或我的笔记“On a Binomial Sampler”。）

以下是二项式 (1/2) 算法的纯 Python 实现，不依赖于底层操作系统的 32 位/64 位支持。

import random
import math

def binomhalf(n):
   if n<4: return sum(random.randint(0,1) for i in range(n))
   if n%2==1: return random.randint(0,1)+binomhalf(n-1)
   m=int(math.sqrt(n))+1
   while True:
      k=0
      while random.randint(0,1)==0: k+=1
      i=k*m+random.randint(0,m-1)
      ret=n//2+i if random.randint(0,1)==0 else n//2-i-1
      if ret<0 or ret>n: continue
      expo=-random.expovariate(1)
      p=math.lgamma(n+1)-math.lgamma(ret+1)-math.lgamma((n-ret)+1)+ \
          math.log(m)+math.log(2)*((k-n)+2)
      if expo<=p: return ret

参考：

K. Bringmann, F. Kuhn 等人，“内部 DLA：物理增长模型的有效模拟”。在：过程中。第 41 届自动机、语言和编程国际学术讨论会 (ICALP'14)，2014 年。
Farach-Colton, M. 和 Tsai, M.T.，2015 年。精确的亚线性二项式采样。算法 73(4)，第 637-651 页。

【讨论】：