支持64位或更高位数的numpy.random.binomial替代方案？

Question

我需要根据二项分布随机生成一系列数字。Numpy的随机数工具可以做到这一点，但不幸的是，它似乎只能处理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在64位平台上确实支持64位，除了Windows；不幸的是，我现在用的是Windows）。

Answer 1

有一种精确且高效的方法可以生成二项分布的随机数，特别是当参数 n 很大时，这种方法在2014年的Bringmann等人的论文中有详细描述。以下是与该论文中描述的算法等效的步骤：

1. 如果 n 小于4，就生成 n 个不偏的随机位（0或1），然后返回它们的和。如果 n 是奇数，则先用 n 减去1来运行这个算法，然后把一个不偏的随机位的值加到结果上，最后返回这个结果。
2. 将 m 设置为 n 的平方根向下取整再加1。
3. （首先，从二项曲线的包络线中采样。）生成不偏的随机位（0或1），直到生成一个0为止。将生成的1的数量记为 k。
4. 随机选择一个在 [0, m) 范围内的整数作为 s，然后将 i 设置为 k*m + s。
5. 将 ret 设置为 n/2 + i 或 n/2 - i - 1，两者的概率相等。
6. （第二步，接受或拒绝 ret。）如果 ret 小于0或大于 n，就回到第3步。
7. 以概率 choose(n, ret)*m*2^{k - n - 2} 返回 ret。否则，回到第3步。（这里的 choose(n, k) 是一个二项系数。下面的Python代码计算这个概率的对数。）

（注意，Bringmann论文中的算法比这更复杂，部分原因是为了避免由于精度有限而导致的舍入误差。此外，Farach-Colton和Tsai在2015年展示了如何将采样二项分布 n, p 的问题简化为采样二项分布 n, 1/2 的问题。想了解更多细节，可以查看这些论文或我的笔记 "关于二项分布采样器"。

以下是一个纯Python实现的二项分布(1/2)算法，不依赖于底层操作系统的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。精确的亚线性二项采样。Algorithmica 73(4), 第637-651页。

Answer 2

在一些64位的机器上，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上遇到问题，可以考虑使用正态近似来处理二项分布。在这么大的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