Re: Poisson distribution

[kanderso%dino.bbn.com@RELAY.CS.NET]

  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)))