A efficient binomial random number generator code

Question

The relevant question is: Algorithm to generate Poisson and binomial random numbers?

I just take her description for the Binomial random number:

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.

There is a trivial solution in Algorithm to generate Poisson and binomial random numbers? through using iterations:

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

However, my purpose of pursing a binomial random number generator is just to avoid the inefficient loops (i from 0 to n). My n could be very large. And p is often very small.

A toy example of my case could be: n=1*10^6, p=1*10^(-7).

The n could range from 1*10^3 to 1*10^10.

Answer 1

If you have small p values, you'll like this one better than the naive implementation you cited. It still loops, but the expected number of iterations is O(np) so it's pretty fast for small p values. If you're working with large p values, replace p with q = 1-p and subtract the return value from n. Clearly, it will be at its worst when p = q = 0.5.

public static int getBinomial(int n, double p) {
   double log_q = Math.log(1.0 - p);
   int x = 0;
   double sum = 0;
   for(;;) {
      sum += Math.log(Math.random()) / (n - x);
      if(sum < log_q) {
         return x;
      }
      x++;
   }
}

The implementation is a variant of Luc Devroye's "Second Waiting Time Method" on page 522 of his text "Non-Uniform Random Variate Generation."

There are faster methods based on acceptance/rejection techniques, but they are substantially more complex to implement.

Answer 2

I could imagine one way to speed it up by a constant factor (e.g. 4).

After 4 throws you will toss a head 0,1,2,3 or 4.

The probabilities for it are something like [0.6561, 0.2916, 0.0486, 0.0036, 0.0001].

Now you can generate one number random number and simulate 4 original throws. If that's not clear how I can elaborate a little more.

This way after some original pre-calculation you can speedup the process almost 4 times. The only requirement for it to be precise is that the granularity of your random generator is at least p^4.