[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

*To*: Gerd Timmermann <GT%feldberg.sger.dialnet.symbolics.com@RELAY.CS.NET>*Subject*: Re: Poisson distribution*From*: kanderso%dino.bbn.com@RELAY.CS.NET*Date*: Wed, 6 Sep 89 14:25:50 EDT*Cc*: slug@RIVERSIDE.SCRC.SYMBOLICS.COM*In-reply-to*: Your message of Tue, 22 Aug 89 19:19:00 +0200. <19890822171913.0.GT@BUNGSBERG.sger.dialnet.symbolics.com>

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 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))) (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**:**Poisson distribution***From:*Gerd Timmermann <GT%feldberg.sger.dialnet.symbolics.com@RELAY.CS.NET>

- Prev by Date:
**Genera 7.2 ECO Tape #1, June 1988** - Next by Date:
**Using TI printer on Symbolics** - Previous by thread:
**Poisson distribution** - Next by thread:
**[no subject]** - Index(es):