A binomial random number generating algorithm that works when $Np$ is very small I need to generate binomial random numbers:


  A binomial random number is the number of heads in $N$ tosses of a coin with probability $p$ of a heads on any single toss. If you generate $N$ uniform random numbers on the interval $(0,1)$ and count the number less than $p$, then the count is a binomial random number with parameters $N$ and $p$.


In my case, my $N$ could range from $10^3$ to $10^{10}$ and my $p$ is around $10^{-7}$.  Often my $Np$ is around $10^{-3}$.

There is a trivial implementation to generate such binomial random number through loops:

getBinomial(int N, double p):
  x = 0
  repeat N times:
    if getUniformRandom() &lt; p: # getUniformRandom() returns a real number in (0,1) 
      x = x+1
  return x


This na&iuml;ve implementation is very slow, $O(N)$. I tried the Acceptance Rejection/Inversion method [1] implemented in the Colt (http://acs.lbl.gov/software/colt/) lib. It is very fast, but the distribution of its generated number only agrees with the na&iuml;ve implementation when $Np$ is not very small. In my case when $Np = 10^{-3}$, the na&iuml;ve implementation can still generate the number 1 after many runs, but the Acceptance Rejection/Inversion method can never generate the number 1 (always returns 0).

Does anyone know what is the problem here? Or can you suggest a better binomial random number generating algorithm that can solve my case?

[1] V. Kachitvichyanukul, B.W. Schmeiser (1988): Binomial random variate generation, Communications of the ACM 31, 216-222.
If $Np$ is small then the distribution is very close to Poisson, which is easier to generate, especially when the expectation is small.

In your specific case, most of the time the result is zero, and the probability that this happens can be calculated exactly. This optimization will probably result in a speed-up of $10^3$.