I found a very simple algorithm that draws values from a Poisson distribution from this project.

The algorithm's code in Java is:

public final int poisson(double a) {
        double limit = Math.exp(-a), prod = nextDouble();
        int n;
        for (n = 0; prod >= limit; n++)
            prod *= nextDouble();
        return n;

nextDouble() is a function from the Random package in Java that returns a uniformly distributed random double, for example 0.885598042879084.

I can't understand how this creates a Poisson distribution.

Can someone explain?

  • $\begingroup$ what is nextDouble()? $\endgroup$
    – PA6OTA
    May 7, 2014 at 16:18
  • $\begingroup$ @PA6OTA Good question, I added the meaning and a link to the API. Sorry for leaving it out initially $\endgroup$ May 7, 2014 at 16:22

4 Answers 4


It is related to the Poisson process: suppose $a$ is fixed and $N$ is the number of independent Exponential (mean 1) RV's to be added until the sum exceeds $a$. In this case, $N \sim Poisson(a)$.

In the code above, everything is anti-logged. For example, instead of adding Exponentials, they multiply Uniforms, due to the relation $Y = - \ln(U)$ is an Exponential (mean 1) RV whenever $U$ is Uniform[0,1].

  • $\begingroup$ doesn't help too much... $\endgroup$ May 7, 2014 at 16:55
  • $\begingroup$ And yet this is exactly the reason why the algorithm works. $\endgroup$
    – Did
    Dec 1, 2014 at 9:28

Suppose that we have $k$ independent and identically distributed exponential random variables $X_1, \dotsc, X_k$ with parameter $\lambda$. If we define a counting process $\{N(t)\}_{t \geq 0}$ such that $S_k := X_1 + \dotsb + X_k$ is the occurring time of the $k$-th event, then this is a Poisson process with rate $\lambda$. This means that $\mathbb{P}\{N(t) \leq k\} = \mathbb{P}\{S_k \geq t\}$, or more informally that you can describe the distribution of $N(t)$ through the distribution of $S_k$.

Now consider that if $X$ is an exponential random variable with parameter $\lambda$, then $Y := F_X(X) = 1 - e^{-\lambda X}$ is a standard uniform random variable (here $F_X$ is the cumulative distribution function of $X$). Note also that $$ F_X^{-1}(Y) = -\frac{1}{\lambda}\log(1 - Y) $$ is an exponential random variable with parameter $\lambda$. Since also $U := 1 - Y$ is a standard uniform random variable, the random variable $$ X' = -\frac{1}{\lambda}\log U $$ is again an exponential random variable with parameter $\lambda$.

We want to count (at most) $k$ events in a given time interval $[0, t]$ (that we can suppose without loss of generality to be $[0, 1]$, i.e. we arbitrarily fix $t = 1$), that is we want $N(1) \leq k$. By what we observed at the beginning, this corresponds to the condition $S_k \geq 1$, which, in its turn, corresponds to the following: $$ \sum_{i = 1}^k X_i \geq 1 \iff \sum_{i = 1}^k -\frac{1}{\lambda}\log U_i \geq 1 \iff -\frac{1}{\lambda}\log\bigg( \prod_{i = 1}^k U_i \bigg) \geq 1 \iff \prod_{i = 1}^k U_i \leq e^{-\lambda}, $$ which is precisely what is described in Knuth's algorithm.


From wikipedia:

The Poisson distribution is a discrete probability distribution that expresses the probability of a given number of events occurring in a fixed interval of time and/or space if these events occur with a known average rate and independently of the time since the last event.

I think that you could get something similar by adding together random values until you hit a maximum, as the random function has a "known average rate" and is basically "independent" of the last function call. It looks to me like multiplying it like so is just a transformation to simplify the code.


Yes, the previous answer is correct. This is simply a version of Knuth's algorithm. It gets slower the larger the parameter lambda.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.