Alternative to numpy.random.binomial that allows 64 bits or more?

Question

I need to randomly generate a series of numbers according to the binomial distribution. Numpy's random suite provides a means to do this, but unfortunately it seems to be limited to handling 32-bit numbers for the value of n and I want to work with values outside that range. 64-bits should suffice, although arbitrarily higher precision would work just as well.

Example output:

>>> 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

Is there an alternative I could use? Or a means to get numpy to use 64-bit numbers internally in this random generator?

Or do I need to saddle up and roll my own?

(As pointed out by Robert Klein, numpy does do 64-bit on 64-bit platforms except Windows; unfortunately I'm using Windows).

Answer 1

On machines where a C long integer is 64-bit, numpy.random.binomial() will accept and generate 64-bit integers. Most 64-bit platforms except for Windows are such. For example, on my 64-bit OS X machine:

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

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

Alternately, if you are stuck on Windows, consider using the Normal Approximation to the binomial distribution. At such large n, the approximation should be excellent.

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