;;; This file contains the procedures for returning a random number
;;; distributed around a mean by a poisson distribution.

(in-package "USER")


(load "random.lisp")


(defvar *lambda-term* 0.0)

;;; This code comes from Melanie Mitchell (via Annie Wu) but I
;;; enhanced it to operate on means > 1.

(defun poisson (mean &optional (compute-lambda-term? t))
   ; Returns a random number (an integer) from a poisson distribution.
   ; The mean of the distribution is passed into the function as well
   ; as an optional flag (compute-lambda-term?).  If the flag is t,
   ; then *lambda-term* is set to expt(- mean).   Otherwise,
   ; *lambda-term* is assumed to have the value exp (- mean).  Use of
   ; this term avoids unnecessary recalculation of the exponent.  Note,
   ; *lambda-term* is not identical with lambda of the poisson
   ; distribution, since lambda = mean.
  
  (let ((uniform (random-unit))
	term
	(k-value 0)
	cum)
    (declare (special *lambda-term*))
    (if (< mean 0) (error "Error in Poisson!!: bad mean"))
  
    (if compute-lambda-term?
	(setf *lambda-term* (exp (- mean))))
  
    ;; Compute the cumulative density function (c.d.f.) for successive
    ;; values of k-value.  When the c.d.f. becomes greater than the
    ;; random variable chosen from the Uniform distribution (uniforms),
    ;; we quit.  Since k must be non-negative, and since
    ;; prob(k=0) = lambda_term, k will start out at 0 with the other values
    ;; initialized accordingly.
  
    (setf term *lambda-term*)
    (setf cum *lambda-term*)
    (loop
     (if (or (> cum uniform)
	     (and (> k-value mean) (<= term 0.0000005))) ; Cutoff for high end
					; of distribution.
	 (return k-value))
  
     (incf k-value)
     (setf term (/ (* term  mean) k-value))
     (setf cum (+ cum term)))
    k-value))
  
  ;------------------------------------------------------------------------

;;; Returns a random fraction over the interaval [0,1).
(defun random-unit () (/ (random most-positive-fixnum) most-positive-fixnum))