Algorithm to generate Poisson and binomial random numbers?

Question

i've been looking around, but i'm not sure how to do it.

i've found this page which, in the last paragraph, says:

A simple generator for random numbers taken from a Poisson distribution is obtained using this simple recipe: if x₁, x₂, ... is a sequence of random numbers with uniform distribution between zero and one, k is the first integer for which the product x₁ · x₂ · ... · x_k+1 < e^-λ

i've found another page describing how to generate binomial numbers, but i think it is using an approximation of poisson generation, which doesn't help me.

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.

i know there are libraries to do it, but i can't use them, only the standard uniform generators provided by the language (java, in this case).

Solution 1

Poisson distribution

Here's how Wikipedia says Knuth says to do it:

init:
     Let L ← e^(−λ), k ← 0 and p ← 1.
do:
     k ← k + 1.
     Generate uniform random number u in [0,1] and let p ← p × u.
while p > L.
return k − 1.

In Java, that would be:

public static int getPoisson(double lambda) {
  double L = Math.exp(-lambda);
  double p = 1.0;
  int k = 0;

  do {
    k++;
    p *= Math.random();
  } while (p > L);

  return k - 1;
}

---

Binomial distribution

Going by chapter 10 of Non-Uniform Random Variate Generation (PDF) by Luc Devroye (which I found linked from the Wikipedia article) gives this:

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

---

Please note

Neither of these algorithms is optimal. The first is O(λ), the second is O(n). Depending on how large these values typically are, and how frequently you need to call the generators, you might need a better algorithm. The paper I link to above has more complicated algorithms that run in constant time, but I'll leave those implementations as an exercise for the reader. :)

Solution 2

For this and other numerical problems the bible is the numerical recipes book.

There's a free version for C here: http://www.nrbook.com/a/bookcpdf.php (plugin required)

Or you can see it on google books: http://books.google.co.uk/books?id=4t-sybVuoqoC&lpg=PP1&ots=5IhMINLhHo&dq=numerical%20recipes%20in%20c&pg=PP1#v=onepage&q=&f=false

The C code should be very easy to transfer to Java.

This book is worth it's weight in gold for lots of numerical problems. On the above site you can also buy the latest version of the book.

Solution 3

Although the answer posted by Kip is perfectly valid for generating Poisson RVs with small rate of arrivals (lambda), the second algorithm posted in Wikipedia Generating Poisson Random variables is better for larger rate of arrivals due to numerical stability.

I faced problems during implementation of one of the projects requiring generation of Poisson RV with very high lambda due to this. So I suggest the other way.

Solution 4

There are several implementations from CERN in the following library (Java code):

http://acs.lbl.gov/~hoschek/colt/

Concerning binomial random numbers, it is based on the paper from 1988 "Binomial Random Variate Generation", that I recommend to you since they use an optimized algorithm.

Regards

Solution 5

you can add this to build.gradle

implementation 'org.kie.modules:org-apache-commons-math:6.5.0.Final'

and use class PoissonDistribution more detail for class PoissonDistribution

View more solutions

Comments

[Akshat Agarwal]

@Kip In your first example, L stood for Lambda, what did P stand for?

[Kip]

@AkshatAgarwal Actually L stands for e^(−lambda). p is just a number between 0 and 1. Upon each iteration, it gets decreased by a random amount (when it is multiplied by a random number between 0 and 1) until you get to the point that p is less than L. The larger lambda is, the smaller (closer to 0) L will be, meaning the more iterations of the loop we will go through on average before p is less than L.

[isapir]

Is there a limit to lambda? When I pass 1000 or 2000 I get an average of 745 in the results. If I pass lower numbers the results seem correct.