A binomial random number generating algorithm that works when $n \times p$ is very small

I need to generate binomial random numbers:

For example, consider 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 $1\times10^{3}$ to $1\times10^{10}$. My p is around $1\times10^{-7}$.

Often my $n\times p$ is around $1 \times 10^{-3}$.

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

public static int getBinomial(int n, double p) {
int x = 0;
for(int i = 0; i < n; i++) {
if(Math.random() < p)
x++;
}
return x;
}

This native implementation is very slow. I tried the Acceptance Rejection/Inversion method  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 native implementation when $n \times p$ is not very small. In my case when $n\times p = 1\times 10^{-3}$, the native 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.

 V. Kachitvichyanukul, B.W. Schmeiser (1988): Binomial random variate generation, Communications of the ACM 31, 216-222.

The following method is theoretically exact, provided we have a "good" random generator urand()$\in[0,1)$. It uses the geometric distribution, i.e. the number of trials until a probability $p$ event occurs for the first time:

int getGeometric(double p) {
return ceil( log( urand() ) / log(1-p) );
}

int getBinomial(int N, double p) {
int count = 0;
for (;;) {
int wait = getGeometric(p);
if (wait > N) return count;
count++;
N -= wait;
}
}

However, if $p$ is really small, you may want to replace log(1-p) with (-p). Also getGeometric exhibits some rounding artefacts when urand()is very small.

The following is roughly based on suggestions in Knuth, Seminumerical Algorithms, 3.4.1.F(2). Many improvements especially for getBeta and getGamma are possible:

#define NOTSOMANY 50
int getBinomial(int N, double p) {
if (N < NOTSOMANY) {
int count = 0;
for (int i=0; i<N; i++) if (urand() < p) count++;
return count;
}
int a = 1+N/2;
int b = N+1-a;
double X = getBeta(a,b);
if (X>p)
return getBinomial(a-1,p/X);
else
return a + getBinomial(b-1,(p-X)(1-X));
}

#define SOMELIMIT  3.0
// must be > 1.0, should be >= 3.0 to make getGamma fast,
// but should also be fairly small to make loop short
double getBeta(double a, double b) {  // a>0, b>0
if (a < SOMELIMIT && b < SOMELIMIT) {
for (;;) {
double Y1 = exp(ln(urand())/a), Y2 =  exp(ln(urand())/b);
if (Y1+Y2 <1) return Y1/(Y1+Y2);
}
} else {
double X1 = getGamma(a), X2 = getGamma(b);
return X1/(X1+X2);
}
}

double getGamma(double a) {  // a > 1
for (;;) {
double Y = tan(Pi*urand());
double X = sqrt(2*a-1)*Y+a-1;
if ( X>0 && urand()
<= (1+Y*Y)*exp((a-1)*ln(X/(a-1))-sqrt(2*a-1)*Y) )
return X;
}
}

Yet another approach: Instead of $N$ trials with $p$, make $N/2$ trials with $1-(1-p)^2=(2-p)p$ to count how many pairs ($M$, say, with $M\approx Np$ if $p$ is small) of consecutive trials have at least one hit. Among these, the pairs with two hits are $(M,p/(2-p))$ binomially distributed

int getBinomial(int N, double p) {
if (N < NOTSOMANY) {
int count = 0;
for (int i=0; i<N; i++) if (urand() < p) count++;
return count;
}
if ( p > 0.5 )
return N-getBinomial(N,1-p);
else {
int M = getBinomial(N/2,p*(2-p));
int result = M + getBinomial(M, p/(2-p));
if (N&1)
if (urand() < p) result++;
return result;
}
}

I didn't investigate that further, but it is definitely useful only if $Np$ is not big (as each step of the result ultimately comes from some count++)

Here's an algorithm that gives a good approximation precisely when $np$ is small. The probability that there are exactly $k$ heads is $\binom nk p^k(1-p)^{n-k}$. Note that $$\frac{\binom n{k+1} p^{k+1}(1-p)^{n-(k+1)}}{\binom nk p^k(1-p)^{n-k}} = \frac p{1-p} \frac{n-k}{k+1},$$ which is small when $np$ is small; hence each term is very large compared to the next term.

Now the conditional probability of having at least $j+1$ heads, given that you have at least $j$ heads, is $$\frac{\sum_{k=j+1}^n \binom nk p^k(1-p)^{n-k}}{\sum_{k=j}^n \binom nk p^k(1-p)^{n-k}} \approx \frac{\binom n{j+1} p^{j+1}(1-p)^{n-(j+1)}}{\binom nj p^j(1-p)^{n-j}} = \frac p{1-p} \frac{n-j}{j+1},$$ since both sums are dominated by their initial terms.

Therefore you can generate the number of heads this way: provisionally assume there are $0$ heads. Then change $0$ to $1$ with probability $\frac p{1-p} n$, else halt. If you get $1$, change $1$ to $2$ with probability $\frac p{1-p} \frac{n-1}2$, else halt. And so on.