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.