Press "Enter" to skip to content

Simulating discrete probability laws


This post is devoted to the classical algorithm for the simulation of a discrete law \( {P=\sum_{a\in A}p_a\delta_a} \) on an at most countable set \( {A} \). One may number the atoms by selecting a bijection \( {\sigma:A\rightarrow\{1,2,\ldots\}} \). Now if we partition the real interval \( {[0,1]} \) into blocs of respective Lebesgue measure \( {p_{\sigma^{-1}(1)},p_{\sigma^{-1}(2)}} \), etc, for instance one may use the intervals

\[ I_1=[0,p_{\sigma^{-1}(1)}[, I_2=[p_{\sigma^{-1}(1)},p_{\sigma^{-1}(1)}+p_{\sigma^{-1}(2)}[,\ldots \]

and if \( {U} \) is a uniform random variable on \( {[0,1]} \), then \( {\mathbb{P}(U\in I_{\sigma(a)})=p_a} \) for every \( {a\in A} \). The classical algorithm goes then as follows: we first generate \( {U} \) using rand, next, if \( {U\leq p_{\sigma^{-1}(1)}} \) then we output \( {\sigma^{-1}(1)} \), elseif \( {U\leq p_{\sigma^{-1}(1)}+p_{\sigma^{-1}(2)}} \), we output \( {\sigma^{-1}(2)} \), etc. The cost of this algorithm is the number of comparisons \( {N} \) needed. This number is random, and follows actually the image law of \( {P} \) by the map \( {\sigma} \). In particular, we have \( {\mathbb{E}(N)=\sum_{a\in A} \sigma(a)p_{a}} \).

Any \( {\sigma} \) minimizing the average cost \( {\mathbb{E}(N)} \) satisfies to \( {p_{\sigma^{-1}(1)}\geq p_{\sigma^{-1}(2)}\geq\cdots} \). The best average cost is obtained by labeling the atoms in decreasing probability order, and its value is the mean of this relabeled discrete law.

For the geometric distribution of arbitrary mean, and for the Poisson distribution of unit mean, the atoms are already labeled decreasingly. For finite discrete distributions, the cost of the reordering is \( {|A|\log|A|} \) with high probability (quick sort, can be done once for all).

One may ask if one can simulate the Poisson distribution of arbitrary mean \( {\lambda} \) from a generator of the Poisson distribution of unit mean. The answer is positive, and relies simply on thinning. Namely, if \( {X_1,\ldots,X_{\lceil\lambda\rceil}} \) are i.i.d. Poisson random variables with unit mean, then their sum \( {S=X_1+\cdots+X_{\lceil\lambda\rceil}} \) is Poisson with mean \( {\lceil\lambda\rceil} \), and conditional on \( {\{S=n\}} \), if \( {B_1,\ldots,B_n} \) are i.i.d. Bernoulli random variables with mean \( {\lambda/\lceil\lambda\rceil} \) then \( {B_1+\cdots+B_S} \) follows the Poisson distribution of mean \( {\lambda} \) (Poisson-Binomial mixtures are Poisson).

I like very much such kind of concrete elementary mathematics, involving basic algorithmic and probabilistic concepts. Try this with your students!


  1. Xi'an 2014-09-26

    But this solution means the computing time grows with λ. I did not check in Devroye (1985) but I seem to remember that there was no uniformly bounded algorithm for simulating Poisson variates of arbitrary mean…

  2. Djalil Chafaï 2014-09-26

    That’s life. Beyond simulation, I like very much the stability of Poisson point processes by thinning. So beautiful…

Leave a Reply

Your email address will not be published. Required fields are marked *