Algorithms: Algorithm 342: generator of random numbers satisfying the Poisson distribution

Abstract

The Poisson distribution gives the probability that px events will occur in a certain interval or volume, where the expected or mean value of events is npx. Applications are described by B. W. Lindgren and G. W. McElrath [1]. For a Monte Carlo calculation we wish to generate numbers px that satisfy the Poisson distribution, that is to find the inverse of the Poisson function. To do this we generate a pseudo-random number in the interval 0, 1 and find the number px such that random _< (probability that the number is pz or less) and random > (the probability that the number is px-1 or less). poisson carlo returns the value-1 to signal that the procedure was cMled with a value of npx < 0 or too large for the precision of the computer. It is the responsibility of the user to test the calculated value if there is any possibility of the occurrence of the error condition. In order to save computing time, values of the Poisson distribution computed at a previous entry for the same value of npx are stored in the own array pson. The previous value of npx is npxl. The actual parameter corresponding to npxl must be a real identifier, not a constant or an expression. Before the first call of poisson carlo the calling program must set npxl to a value # npx. The number of pson elements that were previously computed and stored is computed. If it is desired to save storage space at the expense of computing time, the upper bound 84 of pson may be reduced, but then the limit of computed near the end of the procedure must also be decreased accordingly. The procedure which generates random is preferably algorithm 266 [3] or 294 [2]. It can be called as the actual parameter in the procedure call of poisson carlo. The author thanks Mr. I. D. Hill for numerous suggestions and corrections which greatly improved the algorithm. begin own integer computed; own real pne; own real array pson [0:84]; integer n; real ps; if npx < 0 then go to error; if npx # npxl then begin computed := 0; pne := pson [0] := exp (-npx); if pnc = 0 then go to error; comment pson [0] is the probability that poisson carlo = O. It cannot be zero unless-npx underflows the argument range of procedure exp. For …

Topics

    0 Figures and Tables

      Download Full PDF Version (Non-Commercial Use)