# 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() < p: # getUniformRandom() returns a real number in (0,1)
x = x+1
return x


This naï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ïve implementation when $Np$ is not very small. In my case when $Np = 10^{-3}$, the naï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.

• Is this a question about how to program things in Java, or about algorithms for generating such random numbers (independen of the language)? – Raphael May 10 '14 at 11:08
• @Raphael: I edited to make it more obvious it is an algorithms question. – Wandering Logic May 10 '14 at 13:29
• 1) Floats have finite precision; small numbers can be a problem. 2) Did you excute a goodness-of-fit test on either algorithm? For small $p$, seeing lots of zeroes is not that unlikely. 3) Did you check Knuth's TAoCP? He has a whole chapter on generating random numbers of different distributions. – Raphael May 10 '14 at 14:42
• you should mention and link when you cross-post or repost from other sites. stackoverflow.com/questions/23561551/…, stackoverflow.com/questions/23579630/…, – usul May 11 '14 at 8:39

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