;; -*- Mode:Lisp; Syntax:ANSI-Common-LISP; Coding:us-ascii-unix; fill-column:158 -*-
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;
;; @file      use-prob.lisp
;; @author    Mitch Richling <https://www.mitchr.me>
;; @brief     Augments and supports :mjr_probau.@EOL
;; @std       Common Lisp
;; @see       tst-prob.lisp
;; @copyright
;;  @parblock
;;  Copyright (c) 1997,1998,2004,2010,2011,2012,2015, Mitchell Jay Richling <https://www.mitchr.me> All rights reserved.
;;
;;  Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
;;
;;  1. Redistributions of source code must retain the above copyright notice, this list of conditions, and the following disclaimer.
;;
;;  2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions, and the following disclaimer in the documentation
;;     and/or other materials provided with the distribution.
;;
;;  3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software
;;     without specific prior written permission.
;;
;;  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
;;  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
;;  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
;;  OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
;;  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
;;  DAMAGE.
;;  @endparblock
;; @todo      mjr_prob_bernoulli-pdf: Reconsider "always double-float" rule.  Document why we do it.@EOL@EOL
;; @todo      mjr_prob_negative-binomial-pdf: Optimize!@EOL@EOL
;; @todo      mjr_prob_negative-binomial-pdf: Add a :poisson option for :algorithm.@EOL@EOL
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defpackage :MJR_PROB
  (:USE :COMMON-LISP
        :MJR_NUMU
        :MJR_COMBE
        :MJR_CMP
        :MJR_CHK
        :MJR_PRNG
        :MJR_PROBU)
  (:DOCUMENTATION "Brief: Augments and supports :MJR_PROBAU.;")
  (:EXPORT #:mjr_prob_help
           #:mjr_prob_exponential-pdf        #:mjr_prob_exponential-cdf        #:mjr_prob_exponential-ccdf       #:mjr_prob_exponential-icdf #:mjr_prob_exponential-prng
           #:mjr_prob_std-normal-pdf         #:mjr_prob_std-normal-cdf         #:mjr_prob_std-normal-ccdf        #:mjr_prob_std-normal-icdf  #:mjr_prob_std-normal-prng
           #:mjr_prob_normal-pdf             #:mjr_prob_normal-cdf             #:mjr_prob_normal-ccdf                                        #:mjr_prob_normal-prng
           #:mjr_prob_poisson-pdf            #:mjr_prob_poisson-cdf            #:mjr_prob_poisson-ccdf                                       #:mjr_prob_poisson-prng
           #:mjr_prob_bernoulli-pdf          #:mjr_prob_bernoulli-cdf          #:mjr_prob_bernoulli-ccdf                                     #:mjr_prob_bernoulli-prng
           #:mjr_prob_binomial-pdf           #:mjr_prob_binomial-cdf           #:mjr_prob_binomial-ccdf                                      #:mjr_prob_binomial-prng
           #:mjr_prob_geometric-pdf          #:mjr_prob_geometric-cdf          #:mjr_prob_geometric-ccdf                                     #:mjr_prob_geometric-prng
           #:mjr_prob_negative-binomial-pdf  #:mjr_prob_negative-binomial-cdf  #:mjr_prob_negative-binomial-ccdf                             #:mjr_prob_negative-binomial-prng
           ))

(in-package :MJR_PROB)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_help ()
"Help for MJR_PROB:

See :MJR_PROBU for some vocabulary and notation used in this package.

While this package provides some genuinely useful PDFs; The primary point of this package is to both augment and support :MJR_PROBAU.

It augments it by providing traditionally parametrized versions of some of the PDFs (bernoulli, binomial, geometric, negative-binomial) in :MJR_PROBAU, and it
supports it by providing PDFs (exponential, std-normal, normal, & poisson) useful for approximating the PDFs in :MJR_PROBAU."
  (documentation 'mjr_prob_help 'function))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_exponential-pdf (x mu &key (algorithm :direct))
  "PDF: $\\mu e^{-\\mu x}$"
  (cond ((not (numberp x))               (error "mjr_prob_exponential-pdf: x must be a number!"))
        ((complexp x)                    (error "mjr_prob_exponential-pdf: x must be real (i.e. not complex)!"))
        ((not (numberp mu))              (error "mjr_prob_exponential-pdf: mu must be a number!"))
        ((complexp mu)                   (error "mjr_prob_exponential-pdf: mu must be real (i.e. not complex)!"))
        ((<= mu 0)                       (error "mjr_prob_exponential-pdf: mu must be greater than 0!"))
        ((not (equal algorithm :direct)) (error "mjr_prob_exponential-pdf: Unknown algorithm")))
  (if (< x 0)
      0
      (* mu (exp (* -1 mu x)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_exponential-cdf (x mu &key (algorithm :direct))
  "CDF: $1-e^{-\\mu x}$"
  (cond ((not (numberp x))               (error "mjr_prob_exponential-cdf: x must be a number!"))
        ((complexp x)                    (error "mjr_prob_exponential-cdf: x must be real (i.e. not complex)!"))
        ((not (numberp mu))              (error "mjr_prob_exponential-cdf: mu must be a number!"))
        ((complexp mu)                   (error "mjr_prob_exponential-cdf: mu must be real (i.e. not complex)!"))
        ((<= mu 0)                       (error "mjr_prob_exponential-cdf: mu must be greater than 0!"))
        ((not (equal algorithm :direct)) (error "mjr_prob_exponential-cdf: Unknown algorithm")))
  (- 1 (exp (* -1 mu x))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_exponential-ccdf (x mu &key (algorithm :direct))
  (cond ((not (numberp x))               (error "mjr_prob_exponential-ccdf: x must be a number!"))
        ((complexp x)                    (error "mjr_prob_exponential-ccdf: x must be real (i.e. not complex)!"))
        ((not (numberp mu))              (error "mjr_prob_exponential-ccdf: mu must be a number!"))
        ((complexp mu)                   (error "mjr_prob_exponential-ccdf: mu must be real (i.e. not complex)!"))
        ((<= mu 0)                       (error "mjr_prob_exponential-ccdf: mu must be greater than 0!"))
        ((not (equal algorithm :direct)) (error "mjr_prob_exponential-ccdf: Unknown algorithm")))
  (- 1 (mjr_prob_exponential-cdf x mu)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_exponential-icdf (p mu &key (algorithm :direct))
  "ICDF: $\\frac{1}{\\mu} \\ln\\left(\\frac{-1}{p-1}\\right)$"
  (cond ((not (numberp p))               (error "mjr_prob_exponential-icdf: p must be a number!"))
        ((complexp p)                    (error "mjr_prob_exponential-icdf: p must be real (i.e. not complex)!"))
        ((> p 1)                         (error "mjr_prob_exponential-icdf: p must be less or equal to 1!"))
        ((< p 0)                         (error "mjr_prob_exponential-icdf: p must be greater than or equal to 0!"))
        ((not (numberp mu))              (error "mjr_prob_exponential-icdf: mu must be a number!"))
        ((complexp mu)                   (error "mjr_prob_exponential-icdf: mu must be real (i.e. not complex)!"))
        ((<= mu 0)                       (error "mjr_prob_exponential-icdf: mu must be greater than 0!"))
        ((not (equal algorithm :direct)) (error "mjr_prob_exponential-icdf: Unknown algorithm")))
  (/ (log (/ -1 (- p 1))) mu))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_exponential-prng (mu &key (algorithm :icdf))
  (cond ((not (numberp mu))              (error "mjr_prob_exponential-prng: mu must be a number!"))
        ((complexp mu)                   (error "mjr_prob_exponential-prng: mu must be real (i.e. not complex)!"))
        ((<= mu 0)                       (error "mjr_prob_exponential-prng: mu must be greater than 0!"))
        ((not (equal algorithm :icdf))   (error "mjr_prob_exponential-prng: Unknown algorithm")))
  (mjr_prob_exponential-icdf (mjr_prng_float-oo 0.0 1.0) mu))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_std-normal-pdf (x &key (algorithm :direct))
  "Standard Normal PDF (return is always DOUBLE-FLOAT)

When using IEEE double arithmetic, this function returns zero if x<=-39 or if x>=39.  After x<=-5 or x>=5, the returns are no longer recognized as non-zero by
mjr_chk_!=0 using the default epsilon.

The :ALGORITHM argument must be:
  :DIRECT -- use the classical formula.

Classical formula:
  $$\\frac{e^{\\left(\\frac{x^2}{-2}\\right)}}{\\sqrt{2\\pi}}$$"
  (cond ((not (numberp x))               (error "mjr_prob_std-normal-pdf: X must be a number!"))
        ((complexp x)                    (error "mjr_prob_std-normal-pdf: X must be real (i.e. not complex)!"))
        ((not (equal algorithm :direct)) (error "mjr_prob_std-normal-pdf: Unknown algorithm")))
  (let ((x (float x 1.0d0)))
    (/ (exp (* x x -1/2)) (sqrt (* 2 pi)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_std-normal-cdf (x &key (algorithm :pdf2cdf) (pdf-algorithm :direct))
  (cond ((not (numberp x))                (error "mjr_prob_std-normal-cdf: X must be a number!"))
        ((complexp x)                     (error "mjr_prob_std-normal-cdf: X must be real (i.e. not complex)!"))
        ((not (equal algorithm :pdf2cdf)) (error "mjr_prob_std-normal-cdf: Unknown algorithm")))
  (mjr_probu_pdf2cdf x -6 6 #'mjr_prob_std-normal-pdf nil :algorithm pdf-algorithm))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_std-normal-ccdf (x &key (algorithm :pdf2ccdf) (pdf-algorithm :direct))
  (cond ((not (numberp x))                 (error "mjr_prob_std-normal-ccdf: X must be a number!"))
        ((complexp x)                      (error "mjr_prob_std-normal-ccdf: X must be real (i.e. not complex)!"))
        ((not (equal algorithm :pdf2ccdf)) (error "mjr_prob_std-normal-ccdf: Unknown algorithm")))
  (mjr_probu_pdf2ccdf x -6 6 #'mjr_prob_std-normal-pdf nil :algorithm pdf-algorithm))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_std-normal-icdf (p &key (algorithm :wichura))
  "Inverse CDF

:algorithm is one of
  * :wichura -- Accurate 1 part in 10^16

References:
  Michael Wichura (1988); The Percentage Points of the Normal Distribution; Applied Statistics; Vol 37, Num 3, Pages 477-484; Algorithm AS 241"
  (cond ((not (numberp p))                (error "mjr_prob_std-normal-icdf: P must be a number!"))
        ((complexp p)                     (error "mjr_prob_std-normal-icdf: P must be real (i.e. not complex)!"))
        ((> p 1)                          (error "mjr_prob_std-normal-icdf: P must be less or equal to 1!"))
        ((< p 0)                          (error "mjr_prob_std-normal-icdf: P must be greater than or equal to 0!"))
        ((not (equal algorithm :wichura)) (error "mjr_prob_std-normal-icdf: Unknown algorithm")))
  (let* ((const1 0.180625d+0)
         (const2 1.6d+0)
         (split1 0.425d+00)
         (split2 5.0d+0)
         (q      (- p 0.5d+0)))
    (let ((t1 (- const1 (* q q))))
      (if (<= (abs q) split1)
          (let ((r  t1)
                (a0 3.3871328727963666080d+0)
                (a1 1.3314166789178437745d+2)
                (a2 1.9715909503065514427d+3)
                (a3 1.3731693765509461125d+4)
                (a4 4.5921953931549871457d+4)
                (a5 6.7265770927008700853d+4)
                (a6 3.3430575583588128105d+4)
                (a7 2.5090809287301226727d+3)
                (b0 1.0000000000000000000d+0)
                (b1 4.2313330701600911252d+1)
                (b2 6.8718700749205790830d+2)
                (b3 5.3941960214247511077d+3)
                (b4 2.1213794301586595867d+4)
                (b5 3.9307895800092710610d+4)
                (b6 2.8729085735721942674d+4)
                (b7 5.2264952788528545610d+3))
            (/ (* q (+ (* (+ (* (+ (* (+ (* (+ (* (+ (* (+ (* a7 r) a6) r) a5) r) a4) r) a3) r) a2) r) a1) r) a0))
               (+ (* (+ (* (+ (* (+ (* (+ (* (+ (* (+ (* b7 r) b6) r) b5) r) b4) r) b3) r) b2) r) b1) r) b0)))
          (let ((t2 (if (< q 0) p (- 1 p))))
            (if (<= t2 0)
                (error "ERROR: INFINITY"))
            (let ((t3 (sqrt (- (log t2))))
                  (t4 (if (< q 0) -1.0d0 1.0d0)))
              (if (<= t3 split2)
                  (let ((r  (- t3 const2))
                        (c0 1.42343711074968357734d+0)
                        (c1 4.63033784615654529590d+0)
                        (c2 5.76949722146069140550d+0)
                        (c3 3.64784832476320460504d+0)
                        (c4 1.27045825245236838258d+0)
                        (c5 2.41780725177450611770d-1)
                        (c6 2.27238449892691845833d-2)
                        (c7 7.74545014278341407640d-4)
                        (d0 1.00000000000000000000d+0)
                        (d1 2.05319162663775882187d+0)
                        (d2 1.67638483018380384940d+0)
                        (d3 6.89767334985100004550d-1)
                        (d4 1.48103976427480074590d-1)
                        (d5 1.51986665636164571966d-2)
                        (d6 5.47593808499534494600d-4)
                        (d7 1.05075007164441684324d-9))
                    (* t4 (/ (+ (* (+ (* (+ (* (+ (* (+ (* (+ (* (+ (* c7 r) c6) r) c5) r) c4) r) c3) r) c2) r) c1) r) c0)
                             (+ (* (+ (* (+ (* (+ (* (+ (* (+ (* (+ (* d7 r) d6) r) d5) r) d4) r) d3) r) d2) r) d1) r) d0))))
                  (let ((r   (- t3 split2))
                        (e0  6.65790464350110377720d+0)
                        (e1  5.46378491116411436990d+0)
                        (e2  1.78482653991729133580d+0)
                        (e3  2.96560571828504891230d-1)
                        (e4  2.65321895265761230930d-2)
                        (e5  1.24266094738807843860d-3)
                        (e6  2.71155556874348757815d-5)
                        (e7  2.01033439929228813265d-7)
                        (f0  1.00000000000000000000d+0)
                        (f1  5.99832206555887937690d-1)
                        (f2  1.36929880922735805310d-1)
                        (f3  1.48753612908506148525d-2)
                        (f4  7.86869131145613259100d-4)
                        (f5  1.84631831751005468180d-5)
                        (f6  1.42151175831644588870d-7)
                        (f7  2.04426310338993978564d-15))
                    (* t4 (/ (+ (* (+ (* (+ (* (+ (* (+ (* (+ (* (+ (* e7 r) e6) r) e5) r) e4) r) e3) r) e2) r) e1) r) e0)
                             (+ (* (+ (* (+ (* (+ (* (+ (* (+ (* (+ (* f7 r) f6) r) f5) r) f4) r) f3) r) f2) r) f1) r) f0)))))))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_std-normal-prng (&key (pdf-algorithm :direct) (algorithm :box-muller))
  "Return a pseudo-random number from the standard normal distribution (mean=0 & variance=1)

The :ALGORITHM argument must be:
 * :box-muller    -- 1.0x -- Box-Muller method with cosine. [Box & Muller (1958)] (DEFAULT)
 * :box-mullers   -- 1.0x -- Box-Muller method with sine. [Box & Muller (1958)]
 * :box-mullercs  -- 2.0x -- Box-Muller method with average of cosine & sine [Box & Muller (1958)]
 * :box-mullerp   -- 3.1x -- Box-Muller polar method with averaging.  More stable. [Bell (1968) & Knop (1969)]
 * :box-mullere   -- 2.7x -- Box-Muller exponential method
 * :toms-712      -- 2.6x -- TOMS 712.  Better accuracy, but slower
 * :icdf          -- 4.7x -- Inverse CDF algorithm
 * :accept-reject -- N/A  -- Only returns values in [-6, 6], and is a bit slow
References:
  Box & Muller (1958); A note on the generation of random normal deviates; Annals Math. Stat; Vol 29; pp. 610-611
  Bell (1968); Normal random deviates; Comm. ACM; Algorithm 334
  Knop (1969); Remark on Algorithm 334: normal random deviates; Comm. ACM
  Kinderman & Monahan (1977); ACM Transactions on Mathematical Software; Vol 3 No. 3
  Joseph L. Leva (1992); ACM Transactions on Mathematical Software; Vol 18 No. 4; pp 434-435; TOMS-712"
  (cond ((not (symbolp algorithm)) (error "mjr_prob_std-normal-prng: :ALGORITHM must be a symbol!")))
  (case algorithm
    (:box-muller     (* (sqrt (* -2 (log (mjr_prng_float-co 0.0d0 1.0d0)))) (cos (* 2 pi (random 1.0)))))
    (:icdf           (mjr_prob_std-normal-icdf (mjr_prng_float-oo 0 1)))
    (:toms-712       (loop for urand1 = (mjr_prng_float-co 0.0d0 1.0d0)
                           for urand2 = (* 1.7156 (- (mjr_prng_float-co 0.0d0 1.0d0) 0.5))
                           for x1     = (- urand1 0.449871)
                           for x2     = (- (abs urand2) -0.386595)
                           for qcrv   = (+ (* x1 x1) (* x2 (- (* 0.196 x2) (* 0.25472 x1))))
                           finally (return (/ urand2 urand1))
                           until (< qcrv 0.27597)
                           while (or (> qcrv 0.27846)
                                     (> (* urand2 urand2) (* -4.0 (exp urand1) urand1 urand1)))))
    (:box-mullers    (* (sqrt (* -2 (log (mjr_prng_float-co 0.0d0 1.0d0)))) (sin (* 2 pi (random 1.0)))))
    (:box-mullercs   (let* ((u1 (mjr_prng_float-co 0.0d0 1.0d0))
                            (u2 (mjr_prng_float-co 0.0d0 1.0d0))
                            (t1 (sqrt (* -2 (log u1))))
                            (t2 (* 2 pi u2)))
                       (/ (+ (* t1 (cos t2)) (* t1 (sin t2))) (sqrt 2))))
    (:box-mullerp    (loop for u1   = (mjr_prng_float-cc -1 1)
                           for u2   = (mjr_prng_float-cc -1 1)
                           for u1u2 = (+ (* u1 u1) (* u2 u2))
                           do (if (and (< u1u2 1) (> u1u2 0))
                                  (return (let* ((t1 (sqrt (* -2 (log u1u2) (/ u1u2))))
                                                 (x1 (* t1 u1))
                                                 (x2 (* t1 u2)))
                                            (/ (+ x1 x2) (sqrt 2)))))))
    (:box-mullere    (loop for y1 = (mjr_prob_exponential-prng 1)
                           for y2 = (mjr_prob_exponential-prng 1)
                           do (if (> y2 (/ (expt (- 1 y1) 2) 2))
                                  (return (if (mjr_prng_boolean)
                                              y1
                                              (- y1))))))
    (:accept-reject  (mjr_probu_pdf2prng -6 6 #'mjr_prob_std-normal-pdf nil :algorithm pdf-algorithm))
    (otherwise       (error "mjr_prob_std-normal-prng: Unknown algorithm"))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_normal-pdf (x mean variance &key (algorithm :direct))
  "Normal PDF (return is always DOUBLE-FLOAT)

NOTE: variance = standard deviation squared

Classical formula:
  $$\\frac{e^{\\left(\\frac{(x-\\mu)^2}{-2\\sigma^2}\\right)}}{\\sqrt{2\\pi\\sigma^2}}$$ where $\\mathrm{variance}=\\sigma^2$ and $\\mathrm{mean}=\\mu$"
  (cond ((not (numberp mean))              (error "mjr_prob_normal-pdf: MEAN must be a number!"))
        ((complexp mean)                   (error "mjr_prob_normal-pdf: MEAN must be real (i.e. not complex)!"))
        ((not (numberp variance))          (error "mjr_prob_normal-pdf: VARIANCE must be a number!"))
        ((complexp variance)               (error "mjr_prob_normal-pdf: VARIANCE must be real (i.e. not complex)!"))
        ((not (numberp x))                 (error "mjr_prob_normal-pdf: X must be a number!"))
        ((complexp x)                      (error "mjr_prob_normal-pdf: X must be real (i.e. not complex)!"))
        ((< variance 0)                    (error "mjr_prob_normal-pdf: VARIANCE must be non-negative!"))
        ((not (equal algorithm :direct))   (error "mjr_prob_normal-pdf: Unknown algorithm"))
        ((mjr_cmp_=0 variance 0)           (warn  "mjr_prob_normal-pdf: VARIANCE of zero is silly!")))
  (let ((x         (float x        1.0d0))
        (mean      (float mean     1.0d0))
        (variance  (float variance 1.0d0)))
    (/ (exp (/ (expt (- x mean) 2) (* -2 variance))) (sqrt (* 2 pi variance)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_normal-cdf (x mean variance &key (algorithm :pdf2cdf) (pdf-algorithm :direct))
  (cond ((not (numberp mean))              (error "mjr_prob_normal-cdf: MEAN must be a number!"))
        ((complexp mean)                   (error "mjr_prob_normal-cdf: MEAN must be real (i.e. not complex)!"))
        ((not (numberp variance))          (error "mjr_prob_normal-cdf: VARIANCE must be a number!"))
        ((complexp variance)               (error "mjr_prob_normal-cdf: VARIANCE must be real (i.e. not complex)!"))
        ((not (numberp x))                 (error "mjr_prob_normal-cdf: X must be a number!"))
        ((complexp x)                      (error "mjr_prob_normal-cdf: X must be real (i.e. not complex)!"))
        ((< variance 0)                    (error "mjr_prob_normal-cdf: VARIANCE must be non-negative!"))
        ((not (equal algorithm :pdf2cdf))  (error "mjr_prob_normal-cdf: Unknown algorithm"))
        ((mjr_cmp_=0 variance 0)           (warn  "mjr_prob_normal-cdf: VARIANCE of zero is silly!")))
  (mjr_probu_pdf2cdf x (- mean (* 6 variance)) (+ mean (* 6 variance)) #'mjr_prob_normal-pdf nil mean variance :algorithm pdf-algorithm))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_normal-ccdf (x mean variance &key (algorithm :pdf2ccdf) (pdf-algorithm :direct))
  (cond ((not (numberp mean))              (error "mjr_prob_normal-ccdf: MEAN must be a number!"))
        ((complexp mean)                   (error "mjr_prob_normal-ccdf: MEAN must be real (i.e. not complex)!"))
        ((not (numberp variance))          (error "mjr_prob_normal-ccdf: VARIANCE must be a number!"))
        ((complexp variance)               (error "mjr_prob_normal-ccdf: VARIANCE must be real (i.e. not complex)!"))
        ((not (numberp x))                 (error "mjr_prob_normal-ccdf: X must be a number!"))
        ((complexp x)                      (error "mjr_prob_normal-ccdf: X must be real (i.e. not complex)!"))
        ((< variance 0)                    (error "mjr_prob_normal-ccdf: VARIANCE must be non-negative!"))
        ((not (equal algorithm :pdf2ccdf)) (error "mjr_prob_normal-ccdf: Unknown algorithm"))
        ((mjr_cmp_=0 variance 0)           (warn  "mjr_prob_normal-ccdf: VARIANCE of zero is silly!")))
  (mjr_probu_pdf2ccdf x (- mean (* 6 variance)) (+ mean (* 6 variance)) #'mjr_prob_normal-pdf nil mean variance :algorithm pdf-algorithm))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_normal-prng (mean variance &key (pdf-algorithm :direct) (algorithm :accept-reject) (std-normal-algorithm :box-muller))
  (cond ((not (numberp mean))              (error "mjr_prob_normal-prng: MEAN must be a number!"))
        ((complexp mean)                   (error "mjr_prob_normal-prng: MEAN must be real (i.e. not complex)!"))
        ((not (numberp variance))          (error "mjr_prob_normal-prng: VARIANCE must be a number!"))
        ((complexp variance)               (error "mjr_prob_normal-prng: VARIANCE must be real (i.e. not complex)!"))
        ((< variance 0)                    (error "mjr_prob_normal-prng: VARIANCE must be non-negative!"))
        ((mjr_cmp_=0 variance 0)           (warn  "mjr_prob_normal-prng: VARIANCE of zero is silly!")))
  (case algorithm
    (:accept-reject (mjr_probu_pdf2prng (- mean (* 6 variance)) (+ mean (* 6 variance)) #'mjr_prob_normal-pdf nil mean variance :algorithm pdf-algorithm))
    (:normal        (+ mean (* variance (mjr_prob_std-normal-prng :algorithm std-normal-algorithm))))
    (otherwise      (error "mjr_prob_normal-prng: Unknown algorithm"))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_poisson-pdf (k mu &key (algorithm :direct))
  "Probability (as a DOUBLE-FLOAT) of seeing k events over a time period when the expected number of events over that time is mu.

:algorithm is one of:
   * :direct - Formula may be less stable numerically, but it avoids overflow. approximates the logarithm of the factorial
   * :naive0 - Stable but overflow prone. approximates factorial
   * :naive1 - Stable but overflow prone. exact factorial
   * :normal - Normal approximation good for mu>

Classical formula:
  $$\\frac{\\lambda^k}{k!}\\cdot e^{-\\lambda}$$"
  (cond ((not (numberp mu))        (error "mjr_prob_poisson-pdf: MU must be a number!"))
        ((complexp mu)             (error "mjr_prob_poisson-pdf: MU must be real (i.e. not complex)!"))
        ((< mu 0)                  (error "mjr_prob_poisson-pdf: MU must be non-negative!"))
        ((not (integerp k))        (error "mjr_prob_poisson-pdf: K must be an integer!")))
  (cond ((< k 0)          0)
        ((mjr_chk_!=0 mu) (let ((mu (float mu 1.0d0)))
                            (case algorithm
                              (:direct   (exp (- (* k (log mu)) mu (mjr_numu_log! k))))
                              (:naive0   (/ (* (exp (- mu)) (expt mu k)) (mjr_numu_gamma (1+ k))))
                              (:naive1   (/ (* (exp (- mu)) (expt mu k)) (mjr_combe_! k)))
                              (:normal   (mjr_prob_normal-pdf (+ k 0.5) mu mu))
                              (otherwise (error "mjr_prob_poisson-pdf: Unknown algorithm")))))
        ('t               0))) ;; mu=0 case

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_poisson-cdf (k mu &key (algorithm :pdf2cdf) (pdf-algorithm :direct))
  (cond ((not (numberp mu))                (error "mjr_prob_poisson-cdf: MU must be a number!"))
        ((complexp mu)                     (error "mjr_prob_poisson-cdf: MU must be real (i.e. not complex)!"))
        ((< mu 0)                          (error "mjr_prob_poisson-cdf: MU must be non-negative!"))
        ((not (integerp k))                (error "mjr_prob_poisson-cdf: K must be an integer!"))
        ((not (equal algorithm :pdf2cdf))  (error "mjr_prob_poisson-cdf: Unknown algorithm")))
  (let ((xmax (ceiling (max 300 (cond ((>= mu 100) (+ mu (* 6   mu)))
                                      ((>= mu 50)  (+ mu (* 10  mu)))
                                      ((>= mu 1)   (+ mu (* 200 mu)))
                                      ('t          1))))))
    (mjr_probu_pdf2cdf k 0 xmax #'mjr_prob_poisson-pdf 't mu :algorithm pdf-algorithm)))
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_poisson-ccdf (k mu &key (algorithm :pdf2ccdf) (pdf-algorithm :direct))
  (cond ((not (numberp mu))                (error "mjr_prob_poisson-ccdf: MU must be a number!"))
        ((complexp mu)                     (error "mjr_prob_poisson-ccdf: MU must be real (i.e. not complex)!"))
        ((< mu 0)                          (error "mjr_prob_poisson-ccdf: MU must be non-negative!"))
        ((not (integerp k))                (error "mjr_prob_poisson-ccdf: K must be an integer!"))
        ((not (equal algorithm :pdf2ccdf)) (error "mjr_prob_poisson-ccdf: Unknown algorithm")))
  (let ((xmax (ceiling (max 300 (cond ((>= mu 100) (+ mu (* 6   mu)))
                                      ((>= mu 50)  (+ mu (* 10  mu)))
                                      ((>= mu 1)   (+ mu (* 200 mu)))
                                      ('t          1))))))
    (mjr_probu_pdf2ccdf k 0 xmax #'mjr_prob_poisson-pdf 't mu :algorithm pdf-algorithm)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_poisson-prng (mu &key (algorithm :knuth))
  (cond ((not (numberp mu))             (error "mjr_prob_poisson-prng: MU must be a number!"))
        ((complexp mu)                  (error "mjr_prob_poisson-prng: MU must be real (i.e. not complex)!"))
        ((< mu 0)                       (error "mjr_prob_poisson-prng: MU must be non-negative!"))
        ((not (equal algorithm :knuth)) (error "mjr_prob_poisson-prng: Unknown algorithm")))
  (loop with l = (exp (- mu))
        for k from 1
        for u = (random 1.0)
        for p = 1 then (* p u)
        finally (return (- k 1))
        while (> p l)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_bernoulli-pdf (k p &key (algorithm :direct))
  "Probability of having k successes when trying an experiment 1 time when the probability of success is P

NOTE: If p is a float (single or double), then the result is a DOUBLE-FLOAT

Classical formula:
  $$\\begin{cases}
      1-p & \\text{iff $k=0$} \\\\
      p   & \\text{iff $k=1$} \\\\
      0   & \\text{otherwise} \\\\
    \\end{cases}$$"
  (cond ((not (numberp p))                 (error "mjr_prob_bernoulli-pdf: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_bernoulli-pdf: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_bernoulli-pdf: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_bernoulli-pdf: P must be greater than or equal to 0!"))
        ((not (integerp k))                (error "mjr_prob_bernoulli-pdf: K must be an integer!"))
        ((not (equal algorithm :direct))   (error "mjr_prob_bernoulli-pdf: Unknown algorithm!")))
  (let ((p (mjr_numu_max-accuracy p)))
    (case k
      (0         (- 1 p))
      (1         p)
      (otherwise 0))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_bernoulli-cdf (k p &key (algorithm :pdf2ccdf) (pdf-algorithm :direct))
  (cond ((not (numberp p))                 (error "mjr_prob_bernoulli-cdf: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_bernoulli-cdf: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_bernoulli-cdf: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_bernoulli-cdf: P must be greater than or equal to 0!"))
        ((not (integerp k))                (error "mjr_prob_bernoulli-cdf: K must be an integer!"))
        ((not (equal algorithm :pdf2cdf))  (error "mjr_prob_bernoulli-cdf: Unknown algorithm!")))
  (mjr_probu_pdf2cdf k 0 1 #'mjr_prob_bernoulli-pdf 't p  :algorithm pdf-algorithm))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_bernoulli-ccdf (k p &key (algorithm :pdf2ccdf) (pdf-algorithm :direct))
  (cond ((not (numberp p))                 (error "mjr_prob_bernoulli-ccdf: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_bernoulli-ccdf: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_bernoulli-ccdf: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_bernoulli-ccdf: P must be greater than or equal to 0!"))
        ((not (integerp k))                (error "mjr_prob_bernoulli-ccdf: K must be an integer!"))
        ((not (equal algorithm :pdf2ccdf)) (error "mjr_prob_bernoulli-ccdf: Unknown algorithm!")))
  (mjr_probu_pdf2ccdf k 0 1 #'mjr_prob_bernoulli-pdf 't p  :algorithm pdf-algorithm))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_bernoulli-prng (p &key (algorithm :accept-reject) (pdf-algorithm :direct))
  (cond ((not (numberp p))                 (error "mjr_prob_bernoulli-prng: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_bernoulli-prng: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_bernoulli-prng: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_bernoulli-prng: P must be greater than or equal to 0!")))
  (case algorithm
    (:accept-reject (mjr_probu_pdf2prng 0 1 #'mjr_prob_bernoulli-pdf 't p  :algorithm pdf-algorithm))
    (:bau           (if (< (random 1.0) p) 1 0))
    (otherwise      (error "mjr_prob_bernoulli-prng: Unknown algorithm"))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_geometric-pdf (k p &key (algorithm :direct))
  "The probability of K failures followed by 1 success for K+1 Bernoulli trials (each with success probability of P)

Alternate statement: Number of Bernoulli trial failures required until a success occurs

NOTE: This definition of the geometric distribution is compatible with R in order to facilitate interoperability.

Value of :ALGORITHM determines how the computation is performed.
  * :direct  -- use direct computation using definition

Classical formula:
  $$(1-p)^k\\cdot p$$"
  (cond ((not (numberp p))                 (error "mjr_prob_geometric-pdf: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_geometric-pdf: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_geometric-pdf: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_geometric-pdf: P must be greater than or equal to 0!"))
        ((not (integerp k))                (error "mjr_prob_geometric-pdf: K must be an integer!"))
        ((< k 0)                           (error "mjr_prob_geometric-pdf: K must be non-negative!"))
        ((not (equal algorithm :direct))   (error "mjr_prob_geometric-pdf: Unknown algorithm!")))
  (* p (expt (- 1 p) k)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_geometric-cdf (k p &key (algorithm :pdf2cdf) (pdf-algorithm :direct))
  (cond ((not (numberp p))                 (error "mjr_prob_geometric-cdf: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_geometric-cdf: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_geometric-cdf: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_geometric-cdf: P must be greater than or equal to 0!"))
        ((not (integerp k))                (error "mjr_prob_geometric-cdf: K must be an integer!"))
        ((< k 0)                           (error "mjr_prob_geometric-cdf: K must be non-negative!"))
        ((not (equal algorithm :pdf2cdf))  (error "mjr_prob_geometric-cdf: Unknown algorithm!")))
  (mjr_probu_pdf2cdf k 0 nil #'mjr_prob_geometric-pdf 't p :algorithm pdf-algorithm))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_geometric-ccdf (k p &key (algorithm :pdf2ccdf) (pdf-algorithm :direct))
  (cond ((not (numberp p))                 (error "mjr_prob_geometric-ccdf: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_geometric-ccdf: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_geometric-ccdf: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_geometric-ccdf: P must be greater than or equal to 0!"))
        ((not (integerp k))                (error "mjr_prob_geometric-ccdf: K must be an integer!"))
        ((< k 0)                           (error "mjr_prob_geometric-ccdf: K must be non-negative!"))
        ((not (equal algorithm :pdf2ccdf)) (error "mjr_prob_geometric-ccdf: Unknown algorithm!")))
  (mjr_probu_pdf2ccdf k 0 nil #'mjr_prob_geometric-pdf 't p :algorithm pdf-algorithm))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_geometric-prng (p &key (algorithm :bau))
  (cond ((not (numberp p))                 (error "mjr_prob_geometric-prng: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_geometric-prng: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_geometric-prng: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_geometric-prng: P must be greater than or equal to 0!")))
  (case algorithm
    (:bau           (loop for s = (random 1.0)
                          count (>= s p)
                          until (< s p)))
    (:exponential   (floor (mjr_prob_exponential-prng p)))
    (otherwise      (error "mjr_prob_geometric-pdf: Unknown algorithm"))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_binomial-pdf (k p s &key (algorithm :direct))
  "Probability of exactly K successes in S trials when the probably of success for each trial is P.

Value of :ALGORITHM determines how the computation is performed.
  * :direct  -- use direct computation using definition
  * :poisson -- Use the Poisson approximation
  * :normal  -- Use the Normal approximation (normally good enough
  * :auto    -- Use :poisson when s>100 and p<0.005
                Use :normal when s>1000, p>.2, and p<.8
                Use :direct otherwise

Classical formula:
  $$\\binom{s}{k} p^k (1-p)^{s-k}$$"
  (cond ((not (numberp p))                 (error "mjr_prob_binomial-pdf: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_binomial-pdf: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_binomial-pdf: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_binomial-pdf: P must be greater than or equal to 0!"))
        ((not (integerp s))                (error "mjr_prob_binomial-pdf: S must be an integer!"))
        ((< s 0)                           (error "mjr_prob_binomial-pdf: S must be non-negative!"))
        ((not (integerp k))                (error "mjr_prob_binomial-pdf: K must be an integer!")))
  (let ((algorithm (if (equalp algorithm :auto)
                       (cond ((and (> s 100)  (< p 0.01))        :poisson)
                             ((and (> s 1000) (> p .2) (< p .8)) :normal)
                             ('t                                 :direct))
                       algorithm)))
    (case algorithm
      (:direct   (* (mjr_combe_comb s k) (expt p k) (expt (- 1 p) (- s k))))
      (:poisson  (mjr_prob_poisson-pdf k (* s p)))
      (:normal   (mjr_prob_normal-pdf  k (* s p) (* s p (- 1 p))))
      (otherwise (error "mjr_prob_binomial-pdf: Unknown algorithm")))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_binomial-cdf (k p s &key (algorithm :pdf2cdf) (pdf-algorithm :direct))
  (cond ((not (numberp p))                 (error "mjr_prob_binomial-cdf: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_binomial-cdf: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_binomial-cdf: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_binomial-cdf: P must be greater than or equal to 0!"))
        ((not (integerp s))                (error "mjr_prob_binomial-cdf: S must be an integer!"))
        ((< s 0)                           (error "mjr_prob_binomial-cdf: S must be non-negative!"))
        ((not (equal algorithm :pdf2cdf))  (error "mjr_prob_binomial-cdf: Unknown algorithm!")))
  (mjr_probu_pdf2cdf k 0 s #'mjr_prob_binomial-pdf 't p s :algorithm pdf-algorithm))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_binomial-ccdf (k p s &key (algorithm :pdf2ccdf) (pdf-algorithm :direct))
  (cond ((not (numberp p))                 (error "mjr_prob_binomial-ccdf: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_binomial-ccdf: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_binomial-ccdf: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_binomial-ccdf: P must be greater than or equal to 0!"))
        ((not (integerp s))                (error "mjr_prob_binomial-ccdf: S must be an integer!"))
        ((< s 0)                           (error "mjr_prob_binomial-ccdf: S must be non-negative!"))
        ((not (equal algorithm :pdf2ccdf)) (error "mjr_prob_binomial-ccdf: Unknown algorithm!")))
  (mjr_probu_pdf2ccdf k 0 s #'mjr_prob_binomial-pdf 't p s :algorithm pdf-algorithm))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_binomial-prng (p s &key (algorithm :accept-reject) (pdf-algorithm :direct))
  (cond ((not (numberp p))                 (error "mjr_prob_binomial-prng: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_binomial-prng: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_binomial-prng: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_binomial-prng: P must be greater than or equal to 0!"))
        ((not (integerp s))                (error "mjr_prob_binomial-prng: S must be an integer!"))
        ((< s 0)                           (error "mjr_prob_binomial-prng: S must be non-negative!")))
  (case algorithm
    (:accept-reject (mjr_probu_pdf2prng 0 s #'mjr_prob_binomial-pdf 't p s :algorithm pdf-algorithm))
    (:bau           (loop for i from 1 upto s
                          count (<= (random 1.0) p)))
    (otherwise      (error "mjr_prob_binomial-pdf: Unknown algorithm"))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_negative-binomial-pdf (k p r &key (algorithm :direct))
  "Probability of K failures in a sequence of bernoulli trials with R successes (bernoulli probability is P)

i.e. we do trials till we have R successes, and p(k) is the probability we had k failures in the process.

NOTE: I have used a definition of the negative binomial that is compatible with R in order to facilitate interoperability with R.  It is important to note
      that many sources have K being the successes and R being the failures...

Classical formula:
  $$\\binom{k+r-1}{k}\\cdot (1-p)^k\\cdot p^r$$
Note that:
  $$(-1)^k \\frac{(-r)(-r-1)(-r-2)\\cdots(-r-k+1)}{k!} = (-1)^k\\binom{-r}{k}$$"
  (cond ((not (numberp p))                 (error "mjr_prob_negative-binomial-pdf: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_negative-binomial-pdf: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_negative-binomial-pdf: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_negative-binomial-pdf: P must be greater than or equal to 0!"))
        ((not (integerp r))                (error "mjr_prob_negative-binomial-pdf: R must be an integer!"))
        ((< r 0)                           (error "mjr_prob_negative-binomial-pdf: R must be non-negative!"))
        ((not (integerp k))                (error "mjr_prob_negative-binomial-pdf: K must be an integer!"))
        ((not (equal algorithm :direct))   (error "mjr_prob_negative-binomial-pdf: Unknown algorithm!")))
  (cond ((< k 0) 0)
        ('t      (* (mjr_combe_comb (1- (+ k r)) k)(expt p r) (expt (- 1 p) k)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_negative-binomial-cdf (k p r &key (algorithm :pdf2cdf) (pdf-algorithm :direct))
  (cond ((not (numberp p))                 (error "mjr_prob_negative-binomial-cdf: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_negative-binomial-cdf: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_negative-binomial-cdf: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_negative-binomial-cdf: P must be greater than or equal to 0!"))
        ((not (integerp r))                (error "mjr_prob_negative-binomial-cdf: R must be an integer!"))
        ((< r 0)                           (error "mjr_prob_negative-binomial-cdf: R must be non-negative!"))
        ((not (integerp k))                (error "mjr_prob_negative-binomial-cdf: K must be an integer!"))
        ((not (equal algorithm :pdf2cdf))  (error "mjr_prob_negative-binomial-cdf: Unknown algorithm!")))
  (mjr_probu_pdf2cdf k 0 nil #'mjr_prob_negative-binomial-pdf 't p r :algorithm pdf-algorithm))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_negative-binomial-ccdf (k p r &key (algorithm :pdf2ccdf) (pdf-algorithm :direct))
  (cond ((not (numberp p))                 (error "mjr_prob_negative-binomial-ccdf: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_negative-binomial-ccdf: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_negative-binomial-ccdf: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_negative-binomial-ccdf: P must be greater than or equal to 0!"))
        ((not (integerp r))                (error "mjr_prob_negative-binomial-ccdf: R must be an integer!"))
        ((< r 0)                           (error "mjr_prob_negative-binomial-ccdf: R must be non-negative!"))
        ((not (integerp k))                (error "mjr_prob_negative-binomial-ccdf: K must be an integer!"))
        ((not (equal algorithm :pdf2ccdf)) (error "mjr_prob_negative-binomial-ccdf: Unknown algorithm!")))
  (mjr_probu_pdf2ccdf k 0 nil #'mjr_prob_negative-binomial-pdf 't p r :algorithm pdf-algorithm))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun mjr_prob_negative-binomial-prng (p r &key (algorithm :bau))
  (cond ((not (numberp p))                 (error "mjr_prob_negative-binomial-prng: P must be a number!"))
        ((complexp p)                      (error "mjr_prob_negative-binomial-prng: P must be real (i.e. not complex)!"))
        ((> p 1)                           (error "mjr_prob_negative-binomial-prng: P must be less than or equal to 1!"))
        ((< p 0)                           (error "mjr_prob_negative-binomial-prng: P must be greater than or equal to 0!"))
        ((not (integerp r))                (error "mjr_prob_negative-binomial-prng: R must be an integer!"))
        ((< r 0)                           (error "mjr_prob_negative-binomial-prng: R must be non-negative!"))
        ((not (equal algorithm :bau))      (error "mjr_prob_negative-binomial-prng: Unknown algorithm!")))
  (loop with red  = 0
        with blue = 0
        for i from 1
        finally (return blue)
        do (if (<= (random 1.0) p)
               (incf red)
               (incf blue))
        until (= red r)))