[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Poisson distribution
Date: Tue, 22 Aug 89 19:19+0200
From: Gerd Timmermann <GT@feldberg.sger.dialnet.symbolics.com>
Subject: Poisson distribution
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
the functions you asked for. Please debug them.
;;; 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)))
`(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.
(progn ; Use rejection method.
(multiple-value-bind (sq log-m g)
(let ((log-m (log m)))
(values (sqrt (* 2.0 m)) log-m (- (* m log-m (gammln (+ m 1.0)))))))
;; 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.
(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))))
;;; 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)))
(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)))