# Re: Poisson distribution

```  Date: Tue, 22 Aug 89 19:19+0200
From: Gerd Timmermann <GT@feldberg.sger.dialnet.symbolics.com>
Subject: Poisson distribution
To: slug@riverside.scrc.symbolics.com
Message-ID: <19890822171913.0.GT@BUNGSBERG.sger.dialnet.symbolics.com>

Hello,

I urgently need to generate POISSON distributed random numbers.
Does anyone have a (preferrably LISP) program for that? I know that there are solutions
in IMSL and other commercial FORTRAN packages, but I don't have access to them and
buying would be both to expensive (for this single case) and too time-consuming.

If it isn't too time-consuming, you should look at the book "Numerical
Recipes" by Press et. al.  It is a reasonable source for many
algorithms, and is quite informative.  Here is a LISP translation of

;;;
;;; Numerical Recipes in LISP
;;;

(defmacro memoize (args &body form)
(let ((body `(multiple-value-bind (.entry. .found.) (gethash .args. .table.)
(if .found. (values-list .entry.)
(progn (setf (gethash .args. .table.)
(setq .entry. (multiple-value-list ,@form)))
(values-list .entry.))))))

`(let ((.table. ,(make-hash-table :test (if (cdr args) #'equal #'eq))))
,(if (cdr args)
`(si:with-stack-list (.args. ,@args) ,body)
`(let ((.args. ,(first args))) ,body)))))
;;;
;;; Poisson Deviates p. 207.
;;;

(defun poisson-random (m)
"Returns a integer that is a random deviate drawn from a Poisson
distribution of mean m using RANDOM as thesource of uniform random deviates."
(if (< m 12)
(let ((em -1)
(temp 1.0))		; Multiplying uniform deviates is the
(loop (setf em (1+ em)	; same as adding expoential ones.
temp (* temp (random 1.0)))
(unless (> temp			; Compare exponential rather than
(memoize (m) (exp (- m))))	; computing the log.
(return)))
em)
(progn				; Use rejection method.
(multiple-value-bind (sq log-m g)
(memoize (m)
(let ((log-m (log m)))
(values (sqrt (* 2.0 m)) log-m (- (* m log-m (gammln (+ m 1.0)))))))
(let (y
em)
(loop
;; Y is a deviate from a Lorentzian comparison function.
(setf y  (tan (* pi (random 1.0)))
;; EM is Y shifted and scaled.
em (+ (* sq y) m))
(if (>= em 0.0)			; Keep?
;; The trick for integer-valued distributions.
(progn
(setq em (truncate em))
;; The ratio of the desired distribution to the comparison function.
;; Reject by comparing it to another uniform deviate.
;; The factor 0.9 makes in never exceed 1.0.
(if (<= (random 1.0) (* 0.9 (+ 1.0 (* y y))
(exp (- (* em log-m) (gammln (+ em 1.0)) g))))
(return em))))))))))

;;;
;;; 6.1 Gamma Function, Beta Function, Factorials, Binomial Coefficients.
;;; Page 157.

(defun gammln (x)
"Returns the value of (log (gamma x)) for (> x 0).  Full accuracy is obtained
when (> x 1).  For (< 0 X 1), the reflextion formula (6.1.4) can be used first."
;; Omit double precision if five figure accuracy is good enough.
(let* ((x (- x 1.0d0))
(tmp (let ((tmp (+ x 5.5d0)))
(- (* (+ x 0.5d0) (log tmp)) tmp)))
(ser 1.0d0))
(incf ser (/  76.18009173d0   (incf x 1.0d0)))
(incf ser (/ -86.50532033d0   (incf x 1.0d0)))
(incf ser (/  24.01409822d0   (incf x 1.0d0)))
(incf ser (/  -1.231739516d0  (incf x 1.0d0)))
(incf ser (/   0.120858003d-2 (incf x 1.0d0)))
(incf ser (/  -0.536382d-5    (incf x 1.0d0)))
(+ (log (* 2.50662827465d0 ser)) tmp)))

```

• References: