132 RANDOM NUMBERS 3.4.1 (3) The (Make a web site) Poisson distribution
132 RANDOM NUMBERS 3.4.1 (3) The Poisson distribution with mean p. This distribution is related to the exponential distribution as the binomial distribution is related to the geometric: It represents the number of occurrences, per unit time, of an event that can occur at any instant of time. For example, the number of alpha particles emitted by a radioactive substance in a single second has a Poisson distribution. According to this principle, we can produce a Poisson deviate N by generat- ing independent exponential deviates Xi, Xs, . . . with mean l/p, stopping as soon as Xi + . . . +X, 2 1; then N c m -1. The probability that Xl + . . . + X, > 1 is the probability that a gamma deviate of order m is 2 CL, and this comes to JWm t + eet dt/(m -l)!; hence the probability that N = n is co 1 O tn-le-t & = e-PIli ap tneet dt -n 2 0. s n! (40) If we generate exponential deviates by the logarithm method, the above recipe tells us to stop when -(ln Vi +. . + + In Um)/p 2 1. Simplifying this expression, we see that the desired Poisson deviate can be obtained by calculating e-p, converting it to a fixed point representation, then generating one or more uniform deviates U1, Us, . . . until the product satisfies Vi . . . U, 5 e-p, finally setting Ntm-1. On the average this requires the generation of p + 1 uniform deviates, so it is a very useful approach when fi is not too large. When p is large, we can obtain a method of order log p by using the fact that we know how to handle the gamma and binomial distributions for large orders: First generate X with the gamma distribution of order m = Lap], where Q is a suitable constant. (Since X is equivalent to -ln(Ui . . . Urn), we are essentially bypassing m steps of the previous method.) If X < p, set N c m + Ni, where Ni is a Poisson deviate with mean p -X; and if X > p, set N c Ni, where Ni has the binomial distribution (m -l,p/X). This method is due to J. H. Ahrens and U. Dieter, whose experiments suggest that 6 is a good choice for o. The validity of the above reduction when X 2 p is a consequence of the following important principle: Let Xi, . . . , X, be independent exponential deviates with the same mean; let S, = Xi + . . . + Xj and let Vj = Sj/Sm for 1 5 j 5 m. Then the distribution of VI, VZ, . . . , V,-I is the same as the distribution of m -1 independent uniform deviates sorted into increasing order. To establish this principle formally, we compute the probability that K 5 Ult f I&-I 5 v,-~, given the value of S, = s, for arbitrary values 0 5 Wl 5 .-. 5 21,-l < 1: Let f(vl, ws,. . . , urn-i) be the (m - 1)-fold integral VlS-h/P pevtzlw dtz . . . v,-ls-tl-…-tm–2 J0 pe ,-&Jv2S-tl ,jtmpl ; X J 0 pe-t,-l/~.pe-(S-t’–tm-l)/~ 0