Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

A efficient binomial random number generator code in Java

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.

like image 430
Changwang Zhang Avatar asked May 09 '14 09:05

Changwang Zhang


People also ask

How do you generate a random number in a binomial distribution?

r = binornd( n , p ) generates random numbers from the binomial distribution specified by the number of trials n and the probability of success for each trial p .

How do you generate a random number between 1 and 10 inclusive in Java?

Below is the code showing how to generate a random number between 1 and 10 inclusive. Random random = new Random(); int rand = 0; while (true){ rand = random. nextInt(11); if(rand != 0) break; } System.


2 Answers

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.

like image 149
pjs Avatar answered Sep 22 '22 16:09

pjs


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.

like image 29
Kuba Avatar answered Sep 19 '22 16:09

Kuba