idlastro / Math and Statistics: POIDEV

[Source code]

NAME
POIDEV
PURPOSE
Generate a Poisson random deviate
EXPLANATION
Return an integer random deviate drawn from a Poisson distribution with
a specified mean.    Adapted from procedure of the same name in 
"Numerical Recipes" by Press et al. (1992), Section 7.3
NOTE: This routine became partially obsolete in V5.0 with the 
introduction of the POISSON keyword to the intrinsic functions 
RANDOMU and RANDOMN.     However, POIDEV is still useful for adding 
Poisson noise to an existing image array, for which the coding is much 
simpler than it would be using RANDOMU (see example 1) 
CALLING SEQUENCE
result = POIDEV( xm, [ SEED = ] )
INPUTS
xm - numeric scalar, vector or array, specifying the mean(s) of the 
     Poisson distribution
OUTPUT
result - Long integer scalar or vector, same size as xm
OPTIONAL KEYWORD INPUT-OUTPUT
SEED -  Scalar to be used as the seed for the random distribution.  
        For best results, SEED should be a large (>100) integer.
        If SEED is undefined, then its value is taken from the system 
        clock (see RANDOMU).    The value of SEED is always updated 
        upon output.   This keyword can be used to have POIDEV give 
        identical results on consecutive runs.     
EXAMPLE
(1) Add Poisson noise to an integral image array, im
         IDL> imnoise = POIDEV( im)
(2) Verify the expected mean  and sigma for an input value of 81
         IDL> p = POIDEV( intarr(10000) + 81)   ;Test for 10,000 points
         IDL> print,mean(p),sigma(p)
Mean and sigma of the 10000 points should be close to 81 and 9
METHOD
For small values (< 20) independent exponential deviates are generated 
until their sum exceeds the specified mean, the number of events 
required is returned as the Poisson deviate.   For large (> 20) values,
uniform random variates are compared with a Lorentzian distribution 
function.
NOTES
Negative values in the input array will be returned as zeros.  
REVISION HISTORY
Version 1               Wayne Landsman        July  1992
Added SEED keyword                            September 1992
Call intrinsic LNGAMMA function               November 1994
Converted to IDL V5.0   W. Landsman   September 1997
Use COMPLEMENT keyword to WHERE()        W. Landsman August 2008