(DEFINE-FILE-INFO READTABLE "XCL" PACKAGE (DEFPACKAGE "CDF" ( USE "LISP" "PT")))
(il:filecreated "11-Dec-88 19:25:51" il:{qv}<idl>next>cdf-fns.\;6 25792
il:|changes| il:|to:| (il:vars il:cdf-fnscoms) (il:functions empirical-cdf inverse-empirical-cdf exponential-cdf inverse-exponential-cdf gauss-rv cauchy-rv exponential-rv permutation-rv chisquared-cdf inverse-chisquared-cdf beta-cdf inverse-beta-cdf f-cdf inverse-f-cdf cauchy-cdf inverse-cauchy-cdf logistic-cdf inverse-logistic-cdf test-gamma test-chisquared test-beta gamma log-gamma gamma-cdf inverse-gamma-cdf inverse-gamma-cdf-1 inverse-gamma-cdf-2 gaussian-cdf inverse-gaussian-cdf t-cdf inverse-t-cdf test-cdfs test-t standard-gaussian-cdf horner inverse-standard-gaussian-cdf)
(il:variables *gauss-rv-store* %gamma-coef %inv-gamma-coef-1 %inv-gamma-coef-2 %inv-gamma-coef-3 %inv-gamma-coef-4 least-positive-normalized-single-float %inv-std-gaussian-n1 %inv-std-gaussian-d1 %inv-std-gaussian-n2 %inv-std-gaussian-d2)
(xcl:file-environments "CDF-FNS")
il:|previous| il:|date:| " 8-Dec-88 14:55:49" il:{qv}<idl>next>cdf-fns.\;5)
; Copyright (c) 1988 by Xerox Corporation. All rights reserved.
(il:prettycomprint il:cdf-fnscoms)
(il:rpaqq il:cdf-fnscoms ((il:functions horner) (il:variables least-positive-normalized-single-float) (il:coms (il:* il:|;;| " Gaussian") (il:functions standard-gaussian-cdf) (il:variables %inv-std-gaussian-n1 %inv-std-gaussian-d1 %inv-std-gaussian-n2 %inv-std-gaussian-d2) (il:functions inverse-standard-gaussian-cdf) (il:functions gaussian-cdf inverse-gaussian-cdf)) (il:coms (il:* il:|;;| "T distribution") (il:functions t-cdf inverse-t-cdf)) (il:coms (il:* il:|;;| "Gamma CDF") (il:variables %gamma-coef) (il:functions gamma log-gamma) (il:functions gamma-cdf inverse-gamma-cdf) (il:variables %inv-gamma-coef-1 %inv-gamma-coef-2 %inv-gamma-coef-3 %inv-gamma-coef-4) (il:functions inverse-gamma-cdf-1 inverse-gamma-cdf-2)) (il:coms (il:* il:|;;| "Chi-squared") (il:functions chisquared-cdf inverse-chisquared-cdf)) (il:coms (il:* il:|;;| "Beta") (il:functions beta-cdf inverse-beta-cdf)) (il:coms (il:* il:|;;| "F") (il:functions f-cdf inverse-f-cdf)) (il:coms (il:* il:|;;| "Cauchy") (il:functions cauchy-cdf inverse-cauchy-cdf)) (il:coms (il:* il:|;;| "Logistic ") (il:functions logistic-cdf inverse-logistic-cdf)) (il:coms (il:* il:|;;| "Exponential") (il:functions exponential-cdf inverse-exponential-cdf)) (il:* il:|;;| "Debugging ") (il:functions test-cdfs test-t test-gamma test-chisquared test-beta) (il:* il:|;;| "Random variables") (il:coms (il:* il:|;;| "Gaussian") (il:variables *gauss-rv-store*) (il:functions gauss-rv cauchy-rv exponential-rv)) (il:coms (il:* il:|;;| "Empirical distribution functions") (il:functions empirical-cdf inverse-empirical-cdf) (il:functions permutation-rv)) (eval-when (load) (il:p (export (quote (standard-gaussian-cdf inverse-standard-gaussian-cdf gaussian-cdf inverse-gaussian-cdf t-cdf inverse-t-cdf gamma log-gamma gamma-cdf inverse-gamma-cdf chisquared-cdf inverse-chisquared-cdf beta-cdf inverse-beta-cdf f-cdf inverse-f-cdf cauchy-cdf inverse-cauchy-cdf logistic-cdf inverse-logistic-cdf exponential-cdf inverse-exponential-cdf gauss-rv cauchy-rv exponential-rv empirical-cdf inverse-empirical-cdf permutation-rv)) (find-package "CDF")))) (xcl:file-environments "CDF-FNS")))
(defun horner (x coefficients degree) (il:* il:|;;| "X is a float, COEFFICIENTS is a (pointer) vector of COEFFICIENTS, DEGREE is an integer") (il:* il:|;;| "Degree is one less than the length of COEFFICIENTS.") (do ((result (aref coefficients 0)) (i 1 (1+ i))) ((> i degree) result) (setq result (+ (* result x) (aref coefficients i)))))
(defconstant least-positive-normalized-single-float 1.1754944E-38)
(il:* il:|;;| " Gaussian")
(defun standard-gaussian-cdf (x) (il:* il:|;;;| "Algorithm taken from Appl. Statistics (1973) Vol. 22, No.3") (il:* il:|;;;| "Evaluates the tail area of the standard gaussian density from -infinity to x") (setq x (float x)) (il:* il:|;;| "ltone = (n+9) /3 where n is the number of decimal digits of accuracy. utzero is a value s. that exp{-.5*utzero**2}/ utzero*sqr (2) *pi is just greater than the smallest single-float number") (let ((ltone 5.333334) (utzero 12.9) (fx (abs x)) (change-sign-p (< x 0.0)) value) (declare (type float fx value)) (setq value (if (or (> fx utzero) (and change-sign-p (> fx ltone))) 0.0 (let ((fy (* 0.5 fx fx))) (declare (type float fy)) (if (< fx 1.28) (+ 0.5 (* (- fx) (+ 0.3989423 (/ (* -0.3999035 fy) (+ fy 5.758855 (/ -29.82136 (+ fy 2.624331 (/ 48.696 (+ fy 5.928857))))))))) (/ (* 0.3989423 (exp (- fy))) (+ fx -3.8052E-8 (/ 1.000006 (+ fx 3.980648E-4 (/ 1.986154 (+ fx -0.1516791 (/ 5.293304 (+ fx 4.838592 (/ -15.1509 (+ fx 0.7423809 (/ 30.78993 (+ fx 3.990194)))))))))))))))) (if change-sign-p value (- 1.0 value))))
(defconstant %inv-std-gaussian-n1 (vector -25.44106 41.3912 -18.615 2.506628))
(defconstant %inv-std-gaussian-d1 (vector 3.130829 -21.06224 23.08337 -8.473512 1.0))
(defconstant %inv-std-gaussian-n2 (vector 2.321213 4.850141 -2.297965 -2.78719))
(defconstant %inv-std-gaussian-d2 (vector 1.637068 3.543889 1.0))
(defun inverse-standard-gaussian-cdf (p) (il:* il:|;;;| "Algorithm AS 11, Applied Stat. (1977), Vol. 26, N0. 1") (il:* il:|;;;| "produces standard normal deviate corresponding to lower tail area of p.") (setq p (float p)) (let ((q (- p 0.5))) (declare (type float q)) (if (< (abs q) 0.42) (let ((r (* q q))) (declare (type float r)) (/ (* q (horner r %inv-std-gaussian-n1 3)) (horner r %inv-std-gaussian-d1 4))) (let* ((r (sqrt (- (log (if (> q 0.0) (- 1.0 p) p))))) (value (/ (horner r %inv-std-gaussian-n2 3) (horner r %inv-std-gaussian-d2 2)))) (declare (type float r value)) (if (< q 0.0) (- value) value)))))
(defun gaussian-cdf (x &optional (mu 0.0) (sigma 1.0)) (setq x (float x)) (setq mu (float mu)) (setq sigma (float sigma)) (if (<= sigma 0.0) (error "Sigma must be positive: ~s~%" sigma)) (standard-gaussian-cdf (/ (- x mu) sigma)))
(defun inverse-gaussian-cdf (prob &optional (mu 0.0) (sigma 1.0)) (setq prob (float prob)) (setq mu (float mu)) (setq sigma (float sigma)) (unless (< 0.0 prob 1.0) (error "Probability must lie between 0 and 1: ~s" prob)) (if (<= sigma 0.0) (error "Scale must be positive: ~s" sigma)) (+ mu (* sigma (inverse-standard-gaussian-cdf prob))))
(il:* il:|;;| "T distribution")
(defun t-cdf (x &optional (df 1.0)) (il:* il:|;;;| "adaptation of ACM algorithm 395 (G. W. Hill, v. 13 (1970), 617-619)") (setq x (float x)) (setq df (float df)) (if (<= df 0.0) (error "DF must be greater than zero: ~s" df)) (let* ((zero-tolerance least-positive-normalized-single-float) (1+eps (+ 1.0 single-float-epsilon)) (two-over-pi (/ 2.0 pi)) (t-p x) (t-value (* x x)) (n (truncate df)) (fn df) (z 1.0) (y (/ t-value fn)) (b (+ 1.0 y)) a value) (declare (type float zero-tolerance 1+eps two-over-pi t-p t-value fn z y b a value)) (if (or (>= (abs (- fn (float n))) 1.0E-6) (and (>= n 20) (< t-value fn)) (> n 200)) (progn (il:* il:|;;| "asymptotic series for large or non-integer n") (if (> y 1.0E-6) (setq y (log b))) (setq a (- fn 0.5)) (setq b (* 48.0 a a)) (setq y (* a y)) (setq y (* (+ (/ (+ (/ (+ (* (+ (* (+ (* -0.4 y) -3.3) y) -24.0) y) -85.5) (+ (* 0.8 y y) 100.0 b)) y 3.0) b) 1.0) (sqrt y))) (setq value (* 2.0 (standard-gaussian-cdf (- y))))) (progn (if (and (< n 20) (< t-value 4.0)) (progn (il:* il:|;;| "nested summation of 'cosine' series") (setq a (sqrt y)) (setq y a) (if (eq n 1) (setq a 0.0))) (progn (il:* il:|;;| "'tail' series expansion for large t-values") (setq a (sqrt b)) (setq y (* a fn)) (do ((j 2 (+ j 2))) ((<= (abs (- a z)) zero-tolerance)) (setq z a) (setq y (/ (* y (float (1- j))) (* b (float j)))) (setq a (+ a (/ y (+ fn (float j)))))) (setq fn (+ fn 2.0)) (setq y 0.0) (setq z 0.0) (setq a (- a)))) (do nil ((<= (setq fn (- fn 2.0)) 1+eps)) (setq a (+ y (* a (/ (- fn 1.0) (* b fn)))))) (if (> (abs fn) zero-tolerance) (setq a (* two-over-pi (+ (atan y) (/ a b)))) (setq a (/ a (sqrt b)))) (setq value (- z a)))) (if (>= t-p 0.0) (- 1.0 (* 0.5 value)) (* 0.5 value))))
(defun inverse-t-cdf (prob &optional (df 1.0)) (il:* il:|;;;| "Adapted from the S routine qtdis") (setq prob (float prob)) (setq df (float df)) (if (not (< 0.0 prob 1.0)) (error "Probability must lie between 0 and 1: ~s" prob) (if (< df 1.0) (error "DF must be greater than 1.0: ~s" df))) (let ((half-pi (/ pi 2.0)) (zero-tolerance least-positive-normalized-single-float) (fn df) (pr (* 2.0 (- 1.0 prob))) value) (declare (type float half-pi zero-tolerance fn pr value)) (if (< prob 0.5) (setq pr (* prob 2.0))) (cond ((<= (abs (- fn 2.0)) zero-tolerance) (il:* il:|;;| "For DF = 2") (setq value (sqrt (setq value (- (/ 2.0 (* pr (- 2.0 pr))) 2.0))))) ((<= (abs (- fn 1.0)) zero-tolerance) (il:* il:|;;| "For DF = 1") (setq pr (* pr half-pi)) (setq value (/ (cos pr) (sin pr)))) (t (let* ((a (/ 1.0 (- fn 0.5))) (b (/ 48.0 (* a a))) (c (+ (* (+ (* (+ (/ (* 20700.0 a) b) -98.0) a) -16.0) a) 96.36)) (d (* (+ (/ (+ (/ 94.5 (+ b c)) -3.0) b) 1.0) (sqrt (* a half-pi)) fn)) (x (* d pr)) (y (expt x (/ 2.0 fn)))) (declare (type float a b c d x y)) (if (<= y (+ a 0.05)) (let ((fnp2 (+ fn 2.0))) (declare (type float fnp2)) (setq y (+ (/ (* (+ (* (+ (/ 1.0 (* (+ (/ (+ fn 6.0) (+ fn y)) (* -0.089 d) -0.8220001) fnp2 3.0)) (/ 0.5 (+ fn 4.0))) y) -1.0) (+ fn 1.0)) fnp2) (/ 1.0 y)))) (let ((pr2 (* 0.5 pr))) (declare (type float pr2)) (setq x (inverse-standard-gaussian-cdf pr2)) (setq y (* x x)) (if (< fn 5.0) (setq c (+ c (* 0.3 (- fn 4.5) (+ x 0.6))))) (setq c (+ (* (+ (* (+ (* (+ (* 0.05 d x) -5.0) x) -7.0) x) -2.0) x) b c)) (setq y (* (+ (/ (+ (/ (+ (* (+ (* (+ (* 0.4 y) 6.3) y) 36.0) y) 94.5) c) (- y) -3.0) b) 1.0) x)) (setq y (* a y y)) (if (<= y 0.002) (setq y (+ y (* 0.5 y y))) (setq y (- (exp y) 1.0))))) (setq value (sqrt (* fn y)))))) (if (< prob 0.5) (- value) value)))
(il:* il:|;;| "Gamma CDF")
(defconstant %gamma-coef (vector 0.03586834 -0.1935278 0.4821994 -0.7567041 0.918207 -0.897057 0.988206 -0.5771917 1.0))
(defun gamma (x) (setq x (float x)) (if (<= x 0.0) (error "X must be positive: ~s" x) (if (> x 33.0) (error "X greater than 33: ~s" x))) (let* ((fx (- x (truncate x))) (glfx (if (> fx 0.0) (horner fx %gamma-coef 8) 1.0))) (declare (type float fx glfx)) (cond ((< x 1.0) (/ glfx fx)) ((<= x 2.0) glfx) (t (let ((prod 1.0) (term (+ 1.0 fx)) (test (- x 0.5))) (declare (type float prod term test)) (loop (setq prod (* prod term)) (setq term (+ term 1.0)) (unless (< term test) (return (* prod glfx)))))))))
(defun log-gamma (x) (setq x (float x)) (if (<= x 0.0) (error "X must be positive: ~s" x) (if (<= x 33.0) (log (gamma x)) (let ((x-1 (- x 1.0))) (+ (log (sqrt (* 2.0 pi))) (* (+ x-1 0.5) (log x-1)) (- x-1) (/ 1.0 (* 12.0 x-1)) (/ 1.0 (* 288.0 x-1 x-1)))))))
(defun gamma-cdf (x &optional (eta 1.0)) (setq x (float x)) (setq eta (float eta)) (if (<= x 0.0) (error "X must be positive: ~s" x) (if (<= eta 0.0) (error "Eta must be positive: ~s" eta))) (let ((ceta eta) (term 1.0) (total 1.0)) (declare (type float ceta term total temp)) (if (and (>= x eta) (>= x 7.0)) (do ((i 1 (1+ i)) (limit (truncate (if (< x 11.0) (+ x eta -1.0) (+ eta 10.0))))) ((> i limit) (setq term (* term (- ceta 1.0))) (setq total (+ total (/ term (+ x (- ceta) 2.0)))) (- 1.0 (exp (+ (- x) (* (- eta 1.0) (log x)) (log total) (- (log-gamma eta)))))) (setq ceta (- ceta 1.0)) (setq term (/ (* term ceta) x)) (setq total (+ total term))) (loop (setq ceta (+ ceta 1.0)) (let ((temp (/ x ceta))) (setq term (* term temp)) (setq total (+ total term)) (unless (or (> temp 0.5) (> (/ term total) single-float-epsilon)) (return (/ (exp (setq temp (+ (* eta (log x)) (- x) (log total) (- (log-gamma eta))))) eta))))))))
(defun inverse-gamma-cdf (prob &optional (eta 1.0)) (setq prob (float prob)) (setq eta (float eta)) (if (or (< prob 0.0) (>= prob 1.0)) (error "Prob must lie between 0 and 1: ~s" prob) (if (<= eta 0.0) (error "Eta must be positive: ~s" eta))) (cond ((= prob 0.0) 0.0) ((= eta 0.5) (let ((a (inverse-standard-gaussian-cdf (* 0.5 (- 1.0 prob))))) (* a a 0.5))) ((= eta 1.0) (- (log (- 1.0 prob)))) (t (let ((temp (inverse-gamma-cdf-1 prob eta))) (if (or (>= eta 30.0) (and (>= eta 15.0) (> prob 0.49) (< prob 0.99)) (and (> eta (- 25.2 (* prob 20.8))) (>= prob 0.01) (<= prob 0.99))) temp (inverse-gamma-cdf-2 prob eta temp))))))
(defconstant %inv-gamma-coef-1 (vector -0.01425296 0.01264616))
(defconstant %inv-gamma-coef-2 (vector -0.00588609 0.001400483 -0.01091214 -0.002304527))
(defconstant %inv-gamma-coef-3 (vector -2.728484E-4 0.003135411 -0.009699682 0.01316872 0.002618914 -0.2222222))
(defconstant %inv-gamma-coef-4 (vector 5.406674E-5 3.483789E-5 -7.274761E-4 0.003292181 -0.008729714 0.0 0.4714045 1.0))
(defun inverse-gamma-cdf-1 (prob eta) (if (< eta 0.5) (let* ((a (/ 1.0 (* 9.0 eta))) (b (+ 1.0 (- a) (* (inverse-standard-gaussian-cdf prob) (sqrt a)))) (c (* b b b eta 2.0))) (if (<= c 0.0) (expt prob (* 9.0 a)) c)) (let* ((a (/ 0.5 eta)) (b (* (sqrt a) (il:inv-std-gaussian-cdf prob))) (c (+ (* (+ (* (+ (* (horner b %inv-gamma-coef-1 1) a) (horner b %inv-gamma-coef-2 3)) a) (horner b %inv-gamma-coef-3 5)) a) (horner b %inv-gamma-coef-4 7)))) (* c c c eta))))
(defun inverse-gamma-cdf-2 (prob eta temp) (let* ((gamma-temp (gamma-cdf temp eta)) (p-diff (- prob gamma-temp)) (abs-p-diff (abs p-diff)) (old-abs-correction most-negative-single-float) (old-temp temp) correction abs-correction) (dotimes (i 50) (if (and (>= abs-p-diff gamma-temp) (not (eq i 0))) (setq correction (- (/ correction 2.0))) (progn (setq correction (+ (log abs-p-diff) (- (* (- eta 1.0) (log temp))) temp (log-gamma eta))) (setq correction (if (< p-diff 0.0) (- (exp correction)) (exp correction))) (if (not (eq (< p-diff 0.0) (< correction 0.0))) (setq correction (if (< p-diff 0.0) (- (/ abs-correction 2.0)) (/ abs-correction 2.0)))))) (setq abs-correction (abs correction)) (setq old-abs-correction abs-correction) (setq gamma-temp abs-p-diff) (setq temp (max (+ temp correction) least-positive-normalized-single-float)) (setq old-temp temp) (setq gamma-temp (gamma-cdf temp eta)) (setq p-diff (- prob gamma-temp)) (setq abs-p-diff (abs p-diff)) (unless (and (> abs-p-diff 1.0E-8) (> old-abs-correction (* old-temp 1.0E-6))) (return))) old-temp))
(il:* il:|;;| "Chi-squared")
(defun chisquared-cdf (x &optional (df 1.0)) (setq x (float x)) (setq df (float df)) (cond ((< x 0.0) (error "X must be positive: ~s" x)) ((<= df 0.0) (error "Df must be positive: ~s" df)) ((= x 0.0) 0.0) (t (gamma-cdf (/ x 2.0) (/ df 2.0)))))
(defun inverse-chisquared-cdf (prob &optional (df 1.0)) (setq prob (float prob)) (setq df (float df)) (cond ((not (< 0.0 prob 1.0)) (error "Prob must lie between 0.0 and 1.0: ~s" prob)) ((<= df 0.0) (error "Df must be positive: ~s" df)) ((= prob 0.0) 0.0) (t (* 2.0 (inverse-gamma-cdf prob (/ df 2.0))))))
(il:* il:|;;| "Beta")
(defun beta-cdf (x &optional (a 1.0) (b 1.0)) (il:* il:|;;;| "Source --- This is a transcription of the IMSL routine MDBETA") (setq x (float x)) (setq a (float a)) (setq b (float b)) (cond ((not (< 0.0 x 1.0)) (error "X must lie between 0.0 and 1.0: ~s" x)) ((<= a 0.0) (error "A must be positive: ~s" a)) ((<= b 0.0) (error "B must be positive: ~s" b))) (let ((swap-p (> x 0.5)) temp) (cond ((= x 0.0) 0.0) ((= x 1.0) 1.0) (t (when swap-p (rotatef a b) (setq x (- 1.0 x))) (let ((log-eps (log least-positive-normalized-single-float)) (px (* a (log x))) (pq (log-gamma (+ a b))) (ps (rem b 1.0)) (xint 0.0) (tot 0.0) xb ftemp) (declare (type float single-float-epsilon least-positive-normalized-single-float log-eps px pq ps xint tot xb ftemp)) (if (= ps 0.0) (setq ps 1.0)) (setq xb (- (- (- (+ px (log-gamma (+ ps a))) (log-gamma ps)) (log a)) (log-gamma a))) (il:* il:|;;| "Do Taylor series") (il:* il:|;;| "test for underflow") (when (eq (truncate xb log-eps) 0) (setq xint (* (exp xb) 1.0E+10)) (let ((cnt (* xint a)) (wh 0.0)) (declare (type float cnt wh)) (loop (setq wh (+ wh 1.0)) (setq cnt (* cnt (/ (* (- wh ps) x) wh))) (setq xb (/ cnt (+ a wh))) (setq xint (+ xint xb)) (unless (> (/ xb single-float-epsilon) xint) (return)))) (setq xint (* xint 1.0E-10))) (il:* il:|;;| "Do the integration by parts") (when (> b 1.0) (setq xb (+ px (* b (log (- 1.0 x))) pq (- (log-gamma a)) (- (log b)) (- (log-gamma b)))) (il:* il:|;;| "Scale") (setq ftemp (ftruncate xb log-eps)) (if (< ftemp 0.0) (setq ftemp 0.0)) (setq ps b) (let ((c (/ 1.0 (- 1.0 x))) (cnt (exp (- xb (* ftemp log-eps)))) (wh (- b 1.0))) (declare (type float c cnt wh)) (loop (unless (> wh 0.0) (return)) (setq px (/ (* ps c) (+ a wh))) (if (and (<= px 1.0) (or (<= (/ cnt single-float-epsilon) tot) (<= cnt (/ least-positive-normalized-single-float px)))) (return)) (setq cnt (* cnt px)) (il:* il:|;;| "Rescale") (when (> cnt 1.0) (setq ftemp (- ftemp 1.0)) (setq cnt (* cnt least-positive-normalized-single-float))) (setq ps wh) (if (= ftemp 0.0) (setq tot (+ tot cnt))) (setq wh (- wh 1.0))))) (il:* il:|;;| "Finish up") (setq ftemp (+ tot xint)) (if swap-p (- 1.0 ftemp) ftemp))))))
(defun inverse-beta-cdf (prob &optional (a 1.0) (b 1.0)) (il:* il:|;;;| "Algorithm AS 109, Applied Statist. (1977), Vol. 26, No. 1") (setq prob (float prob)) (setq a (float a)) (setq b (float b)) (cond ((not (< 0.0 prob 1.0)) (error "Prob must lie between 0.0 and 1.0: ~s" prob)) ((<= a 0.0) (error "A must be positive: ~s" a)) ((<= b 0.0) (error "B must be positive: ~s" b))) (if (or (= prob 0.0) (= prob 1.0) (and (= a 1.0) (= b 1.0))) prob (let ((swap-p (> prob 0.5)) (beta (- (+ (log-gamma a) (log-gamma b)) (log-gamma (+ a b)))) value) (declare (type beta value)) (il:* il:|;;| "change tail if necessary") (when swap-p (rotatef a b) (setq prob (- 1.0 prob))) (il:* il:|;;| "Calculate the initial approximation") (let ((r (sqrt (- (log (* prob prob))))) y t1 t2 s h w) (declare (type float r y t1 t2 s w)) (setq y (- r (/ (+ 2.30753 (* 0.27061 r)) (+ 1.0 (* (+ 0.99229 (* 0.04481 r)) r))))) (if (or (<= a 1.0) (<= b 1.0)) (progn (setq r (+ b b)) (setq t1 (/ 1.0 (* 9.0 b))) (setq t2 (+ 1.0 (- t1) (* y (sqrt t1)))) (setq t1 (* r t2 t2 t2)) (if (< t1 0.0) (setq value (- 1.0 (exp (/ (log (+ (* (- 1.0 prob) b) beta)) b)))) (progn (setq t1 (/ (+ (* 4.0 a) r -2.0) t1)) (if (< t1 1.0) (setq value (exp (/ (+ (log (* prob a)) beta) a))) (setq value (- 1.0 (/ 2.0 (+ t1 1.0)))))))) (progn (setq r (/ (+ (* y y) -3.0) 6.0)) (setq s (/ 1.0 (+ a a -1.0))) (setq t1 (/ 1.0 (+ b b -1.0))) (setq h (/ 2.0 (+ s t1))) (setq w (- (/ (* y (sqrt (+ h r))) h) (* (- t1 s) (+ r (/ 5.0 6.0) (/ -2.0 (* 3.0 h)))))) (setq value (/ a (+ a (* b (exp (+ w w))))))))) (il:* il:|;;| "Solve for X by a modified newton-raphson method") (cond ((< value 1.0E-4) (setq value 1.0E-4)) ((> value 0.9999) (setq value 0.9999))) (let ((r (- 1.0 a)) (yprev 0.0) (sq 1.0) (prev 1.0) y oldvalue tx) (declare (type float r yprev sq prev y oldvalue tx)) (loop (setq oldvalue value) (setq y (* (- (beta-cdf value a b) prob) (exp (+ beta (* r (log value)) (* (- 1.0 b) (log (- 1.0 value))))))) (if (<= (* y yprev) 0.0) (setq prev sq)) (let ((g 3.0) adj) (declare (type float g adj)) (loop (loop (setq g (/ g 3.0)) (setq adj (* g y)) (setq sq (* adj adj)) (setq tx (- value adj)) (unless (or (>= sq prev) (< tx 0.0) (> tx 1.0)) (return))) (unless (and (>= prev single-float-epsilon) (>= (* y y) prev) (or (= tx 0.0) (= tx 1.0))) (return)))) (setq value tx) (setq yprev y) (unless (and (>= prev single-float-epsilon) (>= (* y y) single-float-epsilon) (not (= tx oldvalue))) (return)))) (if swap-p (- 1.0 value) value))))
(il:* il:|;;| "F")
(defun f-cdf (x &optional (p 1.0) (q 1.0)) (setq x (float x)) (setq p (float p)) (setq q (float q)) (cond ((< x 0.0) (error "X must be positive: ~s" x)) ((<= p 0.0) (error "P must be positive: ~s" p)) ((<= q 0.0) (error "Q must be positive: ~s" q)) ((= x 0.0) 0.0) (t (let ((temp (* x (/ p q)))) (beta-cdf (/ temp (+ 1.0 temp)) (/ p 2.0) (/ q 2.0))))))
(defun inverse-f-cdf (prob &optional (p 1.0) (q 1.0)) (setq prob (float prob)) (setq p (float p)) (setq q (float q)) (cond ((not (< 0.0 prob 1.0)) (error "Prob must lie between 0.0 and 1.0: ~s" prob)) ((<= p 0.0) (error "P must be positive: ~s" p)) ((<= q 0.0) (error "Q must be positive: ~s" q)) ((= prob 0.0) 0.0) (t (let ((temp (inverse-beta-cdf prob (/ p 2.0) (/ q 2.0)))) (* (/ q p) (/ temp (- 1.0 temp)))))))
(il:* il:|;;| "Cauchy")
(defun cauchy-cdf (x &optional (mu 0.0) (sigma 1.0)) (setq x (float x)) (setq mu (float mu)) (setq sigma (float sigma)) (if (<= sigma 0.0) (error "Sigma must be positive: ~s" sigma) (+ (/ (atan (- x mu) sigma) pi) 0.5)))
(defun inverse-cauchy-cdf (prob &optional (mu 0.0) (sigma 1.0)) (setq prob (float prob)) (setq mu (float mu)) (setq sigma (float sigma)) (cond ((not (< 0.0 prob 1.0)) (error "Prob must lie between 0.0 and 1.0: ~s" prob)) ((<= sigma 0.0) (error "Sigma must be positive: ~s" sigma)) (t (+ mu (* sigma (tan (* pi (+ prob -0.5))))))))
(il:* il:|;;| "Logistic ")
(defun logistic-cdf (x &optional (mu 0.0) (sigma 1.0)) (setq x (float x)) (setq mu (float mu)) (setq sigma (float sigma)) (if (<= sigma 0.0) (error "Sigma must be positive: ~s" sigma) (/ 1.0 (+ 1.0 (exp (/ (- mu x) sigma))))))
(defun inverse-logistic-cdf (prob &optional (mu 0.0) (sigma 1.0)) (setq prob (float prob)) (setq mu (float mu)) (setq sigma (float sigma)) (cond ((not (< 0.0 prob 1.0)) (error "Prob must lie between 0.0 and 1.0: ~s" prob)) ((<= sigma 0.0) (error "Sigma must be positive: ~s" sigma)) (t (+ mu (* sigma (log (/ prob (- 1.0 prob))))))))
(il:* il:|;;| "Exponential")
(defun exponential-cdf (x &optional (mu 1.0)) (setq x (float x)) (setq mu (float mu)) (cond ((< x 0.0) (error "X must be positive: ~s" x)) ((<= mu 0.0) (error "Mu must be positive: ~s" mu)) (t (- 1.0 (exp (- (/ x mu)))))))
(defun inverse-exponential-cdf (prob &optional (mu 1.0)) (setq prob (float prob)) (setq mu (float mu)) (cond ((not (< 0.0 prob 1.0)) (error "Prob must lie between 0.0 and 1.0: ~s" prob)) ((<= mu 0.0) (error "Mu must be positive: ~s" mu)) (t (- (* mu (log (- 1.0 prob)))))))
(il:* il:|;;| "Debugging ")
(defun test-cdfs (cdf inv-cdf &optional (n 100)) (let ((inc (/ 1.0 (1+ n)))) (do ((x inc (+ x inc)) (i 0 (1+ i))) ((eq i n)) (let* ((v (funcall inv-cdf x)) (y (funcall cdf v))) (format t "x ~s inv-cdf ~s cdf ~s diff ~s~%" x v y (abs (- x y)))))))
(defun test-t (&optional (min 1.0) (max 100.0) (inc 2.0)) (do ((df min (+ df inc))) ((> df max)) (format t "df is ~s~%" df) (test-cdfs (function (lambda (x) (t-cdf x df))) (function (lambda (x) (inverse-t-cdf x df))) 10)))
(defun test-gamma (&optional (min 0.2) (max 10.0) (inc 0.2)) (do ((eta min (+ eta inc))) ((> eta max)) (format t "Eta is ~s~%" eta) (test-cdfs (function (lambda (x) (gamma-cdf x eta))) (function (lambda (x) (inverse-gamma-cdf x eta))) 10)))
(defun test-chisquared (&optional (min 1.0) (max 100.0) (inc 3.0)) (do ((df min (+ df inc))) ((> df max)) (format t "Df is ~s~%" df) (test-cdfs (function (lambda (x) (chisquared-cdf x df))) (function (lambda (x) (inverse-chisquared-cdf x df))) 10)))
(defun test-beta (&key (min-a 1.0) (max-a 20.0) (inc-a 2.0) (min-b 1.0) (max-b 20.0) (inc-b 2.0) (n 10)) (do ((a min-a (+ a inc-a))) ((> a max-a)) (format t "A is ~s~%" a) (do ((b min-b (+ b inc-b))) ((> b max-b)) (format t "B is ~s~%" b) (let ((inc (/ 1.0 (1+ n)))) (do ((x inc (+ x inc)) (i 0 (1+ i))) ((eq i n)) (let* ((v (inverse-beta-cdf x a b)) (y (beta-cdf v a b))) (format t "x ~s inv-cdf ~s cdf ~s diff ~s~%" x v y (abs (- x y)))))))))
(il:* il:|;;| "Random variables")
(il:* il:|;;| "Gaussian")
(defvar *gauss-rv-store* nil)
(defun gauss-rv (&optional (mu 0.0) (sigma 1.0)) (il:* il:|;;| "Gaussian variates ala BOX-MUELLER") (setq mu (float mu)) (setq sigma (float sigma)) (if (<= sigma 0.0) (error "SIGMA must be positive: ~s" sigma)) (+ mu (* sigma (let ((store *gauss-rv-store*)) (if store (progn (setq *gauss-rv-store* nil) store) (let ((twopi (* 2.0 pi)) (r1 (random 1.0)) (r2 (random 1.0)) a b) (setq a (sqrt (* -2.0 (log r1)))) (setq b (* twopi r2)) (setq *gauss-rv-store* (* a (cos b))) (* a (sin b))))))))
(defun cauchy-rv (&optional (mu 0.0) (sigma 1.0)) (inverse-cauchy-cdf (random 1.0) mu sigma))
(defun exponential-rv (&optional (mu 0.0)) (inverse-exponential-cdf (random 1.0) mu))
(il:* il:|;;| "Empirical distribution functions")
(defun empirical-cdf (x sorted-vector predicate) (let* ((length (length sorted-vector)) (index (sort:binary-search x sorted-vector predicate :interpolate-p t))) (/ (1+ index) (float length))))
(defun inverse-empirical-cdf (prob sorted-vector) (il:* il:|;;| "Least number such that P( X <= x) >= Prob") (if (not (< 0.0 prob 1.0)) (error "Prob must be between 0.0 and 1.0: ~s" prob)) (aref sorted-vector (1- (ceiling (* prob (length sorted-vector))))))
(defun permutation-rv (vector) (il:* il:|;;;| "Generate a random permutation of VECTOR") (do ((new-vector (copy-seq vector)) (j (1- (length vector)) (1- j))) ((< j 0) new-vector) (let ((i (random (1+ j)))) (rotatef (aref new-vector i) (aref new-vector j)))))
(eval-when (load)
(export (quote (standard-gaussian-cdf inverse-standard-gaussian-cdf gaussian-cdf inverse-gaussian-cdf t-cdf inverse-t-cdf gamma log-gamma gamma-cdf inverse-gamma-cdf chisquared-cdf inverse-chisquared-cdf beta-cdf inverse-beta-cdf f-cdf inverse-f-cdf cauchy-cdf inverse-cauchy-cdf logistic-cdf inverse-logistic-cdf exponential-cdf inverse-exponential-cdf gauss-rv cauchy-rv exponential-rv empirical-cdf inverse-empirical-cdf permutation-rv)) (find-package "CDF"))
)
(xcl:define-file-environment "CDF-FNS" :package (xcl:defpackage "CDF" (:use "LISP" "PT")) :readtable "XCL" :compiler :compile-file)
(il:putprops il:cdf-fns il:copyright ("Xerox Corporation" 1988))
(il:declare\: il:dontcopy
(il:filemap (nil)))
il:stop