(DEFINE-FILE-INFO READTABLE "INTERLISP" PACKAGE "INTERLISP")
(FILECREATED " 7-Dec-88 11:36:26" {QV}<IDL>KOTO>CDFS.;3 26819
changes to%: (VARS CDFSCOMS) (MACROS HORNER)
previous date%: "23-Jun-86 14:01:48" {QV}<IDL>KOTO>CDFS.;1)
(* "
Copyright (c) 1986, 1988 by Xerox Corporation. All rights reserved.
")
(PRETTYCOMPRINT CDFSCOMS)
(RPAQQ CDFSCOMS ((DECLARE%: DONTEVAL@LOAD DOEVAL@COMPILE DONTCOPY (FILES UNBOXEDOPS)) (* ;;; "File created by Coms Manager.") (FNS BETA-CDF CAUCHY-CDF CHISQUARED-CDF EXPONENTIAL-CDF F-CDF GAMLOG GAMMA GAMMA-CDF GAUSSIAN-CDF INV-BETA-CDF INV-CAUCHY-CDF INV-CHISQUARED-CDF INV-EXPONENTIAL-CDF INV-F-CDF INV-GAMMA-CDF INV-GAMMA-CDF1 INV-GAMMA-CDF2 INV-GAUSSIAN-CDF INV-LOGISTIC-CDF INV-LOGNORMAL-CDF INV-STD-GAUSSIAN-CDF INV-T-CDF INV-UNIFORM-CDF LOGISTIC-CDF LOGNORMAL-CDF STD-GAUSSIAN-CDF T-CDF UNIFORM-CDF) (MACROS HORNER) (CONSTANTS (\IEEE.SHORT.MACHEPS (EXPT 2.0 -23)) (\IEEE.SHORT.LEAST.NORMAL (EXPT 2.0 -126)) (\IEEE.SHORT.BASE 2) (\IEEE.SHORT.MANTISSA 23) (\IEEE.SHORT.EXPONENT 8) (\IEEE.SHORT.LEAST.EXPONENT -126) (\IEEE.SHORT.GREATEST.EXPONENT 127) (\IEEE.SHORT.EXPONENT.BIAS 127) (\IEEE.SHORT.PI 3.141593)) (DECLARE%: DONTEVAL@LOAD DOEVAL@COMPILE DONTCOPY (LOCALVARS . T))))
(DECLARE%: DONTEVAL@LOAD DOEVAL@COMPILE DONTCOPY
(FILESLOAD UNBOXEDOPS)
)
(* ;;; "File created by Coms Manager.")
(DEFINEQ
(BETA-CDF
(LAMBDA (X A B) (* jop%: "23-Jun-86 13:30") (* * Source - This is a transcription of the IMSL routine MDBETA) (SETQ X (FLOAT X)) (SETQ A (FLOAT (OR A 1.0))) (SETQ B (FLOAT (OR B 1.0))) (if (OR (LESSP X 0.0) (GREATERP X 1.0)) then (ERROR "X must lie between 0 and 1." X) elseif (LEQ A 0.0) then (ERROR "A must be positive." A) elseif (LEQ B 0.0) then (ERROR "B must be positive." B)) (LET ((SWAP-P (FGREATERP X 0.5)) TEMP) (if (FEQP X 0.0) then 0.0 elseif (FEQP X 1.0) then 1.0 else (if SWAP-P then (* Swap args) (LET ((AA A)) (SETQ A B) (SETQ B AA) (SETQ X (FDIFFERENCE 1.0 X)))) (LET ((EPS (CONSTANT \IEEE.SHORT.MACHEPS)) (EPS1 (CONSTANT \IEEE.SHORT.LEAST.NORMAL)) (ALEPS (CONSTANT (LOG \IEEE.SHORT.LEAST.NORMAL))) (PX (FTIMES A (LOG X))) (PQ (GAMLOG (FPLUS A B))) (PS (FREMAINDER B 1.0)) (XINT 0.0) (TOT 0.0) XB FTEMP) (DECLARE (TYPE FLOATP EPS EPS1 ALEPS PX PQ PS XINT TOT XB FTEMP)) (if (UFEQP PS 0.0) then (SETQ PS 1.0)) (SETQ XB (FDIFFERENCE (FDIFFERENCE (FDIFFERENCE (FPLUS PX (GAMLOG (FPLUS PS A))) (GAMLOG PS)) (LOG A)) (GAMLOG A))) (* * Do Taylor series) (* test for underflow) (if (EQ (UFIX (FQUOTIENT XB ALEPS)) 0) then (SETQ XINT (FTIMES (ANTILOG XB) 1.0E+10)) (bind (CNT ← (FTIMES XINT A)) (WH ← 0.0) declare (TYPE FLOATP CNT WH) repeatwhile (UFGREATERP (FQUOTIENT XB EPS) XINT) do (SETQ WH (FPLUS WH 1.0)) (SETQ CNT (FTIMES CNT (FQUOTIENT (FTIMES (FDIFFERENCE WH PS) X) WH))) (SETQ XB (FQUOTIENT CNT (FPLUS A WH))) (SETQ XINT (FPLUS XINT XB))) (SETQ XINT (FTIMES XINT 1.0E-10))) (* * Do the integration by parts) (if (FGREATERP B 1.0) then (SETQ XB (FPLUS PX (FTIMES B (LOG (FDIFFERENCE 1.0 X))) PQ (FMINUS (GAMLOG A)) (FMINUS (LOG B)) (FMINUS (GAMLOG B)))) (* Scale) (SETQ FTEMP (FLOAT (UFIX (FQUOTIENT XB ALEPS)))) (if (UFLESSP FTEMP 0.0) then (SETQ FTEMP 0.0)) (SETQ PS B) (bind (C ← (FQUOTIENT 1.0 (FDIFFERENCE 1.0 X))) (CNT ← (ANTILOG (FDIFFERENCE XB (FTIMES FTEMP ALEPS)))) (WH ← (FDIFFERENCE B 1.0)) declare (TYPE FLOATP C CNT WH) while (UFGREATERP WH 0.0) do (SETQ PX (FQUOTIENT (FTIMES PS C) (FPLUS A WH))) (if (AND (UFLEQ PX 1.0) (OR (UFLEQ (FQUOTIENT CNT EPS) TOT) (UFLEQ CNT (FQUOTIENT EPS1 PX)))) then (RETURN)) (SETQ CNT (FTIMES CNT PX)) (* Rescale) (if (UFGREATERP CNT 1.0) then (SETQ FTEMP (FDIFFERENCE FTEMP 1.0)) (SETQ CNT (FTIMES CNT EPS1))) (SETQ PS WH) (if (UFEQP FTEMP 0.0) then (SETQ TOT (FPLUS TOT CNT))) (SETQ WH (FDIFFERENCE WH 1.0)))) (* * Finish up) (SETQ FTEMP (FPLUS TOT XINT)) (if SWAP-P then (FDIFFERENCE 1.0 FTEMP) else FTEMP)))))
)
(CAUCHY-CDF
(LAMBDA (X MU SIGMA) (* jop%: "30-May-86 13:12") (* *) (SETQ X (FLOAT X)) (SETQ MU (FLOAT (OR MU 0.0))) (SETQ SIGMA (FLOAT (OR SIGMA 1.0))) (if (LEQ SIGMA 0.0) then (ERROR "Sigma must be positive." SIGMA) else (FPLUS (FQUOTIENT (ARCTAN2 (FDIFFERENCE X MU) SIGMA T) \IEEE.SHORT.PI) 0.5)))
)
(CHISQUARED-CDF
(LAMBDA (X DF) (* jop%: "23-Jun-86 13:30") (* *) (SETQ X (FLOAT X)) (SETQ DF (FLOAT (OR DF 1.0))) (if (FLESSP X 0.0) then (ERROR "X must be positive." X) elseif (LEQ DF 0.0) then (ERROR "DegreesOfFreedom must be positive." DF) elseif (FEQP X 0.0) then 0.0 else (GAMMA-CDF (QUOTIENT X 2.0) (QUOTIENT DF 2.0))))
)
(EXPONENTIAL-CDF
(LAMBDA (X MU) (* jop%: "30-May-86 13:15") (* *) (SETQ X (FLOAT X)) (SETQ MU (FLOAT (OR MU 1.0))) (if (FLESSP X 0.0) then (ERROR "X must be positive" X) elseif (LEQ MU 0.0) then (ERROR "Mu must be positive." MU) else (FDIFFERENCE 1.0 (ANTILOG (FMINUS (FQUOTIENT X MU))))))
)
(F-CDF
(LAMBDA (X P Q) (* jop%: " 2-Jun-86 20:40") (* *) (SETQ X (FLOAT X)) (SETQ P (FLOAT (OR P 1.0))) (SETQ Q (FLOAT (OR Q 1.0))) (if (FLESSP X 0.0) then (ERROR "X must be positive." X) elseif (LEQ P 0.0) then (ERROR "NumeratorDF must be positive." P) elseif (LEQ Q 0.0) then (ERROR "DenominatorDF must be positive." Q) elseif (FEQP X 0.0) then 0.0 elseif (LET ((TEMP (TIMES X (QUOTIENT P Q)))) (BETA-CDF (QUOTIENT TEMP (PLUS 1.0 TEMP)) (QUOTIENT P 2.0) (QUOTIENT Q 2.0)))))
)
(GAMLOG
(LAMBDA (X) (* jop%: "10-Jun-86 14:32") (* *) (SETQ X (FLOAT X)) (if (LEQ X 0.0) then (ERROR "X must be positive." X) elseif (LEQ X 33.0) then (LOG (GAMMA X)) else (LET ((XLESSONE (FDIFFERENCE X 1.0))) (FPLUS (CONSTANT (LOG (SQRT (FTIMES 2.0 \IEEE.SHORT.PI)))) (FTIMES (FPLUS XLESSONE 0.5) (LOG XLESSONE)) (FMINUS XLESSONE) (FQUOTIENT 1.0 (FTIMES 12.0 XLESSONE)) (FQUOTIENT 1.0 (FTIMES 288.0 XLESSONE XLESSONE))))))
)
(GAMMA
(LAMBDA (X) (* jop%: "10-Jun-86 14:29") (* *) (SETQ X (FLOAT X)) (if (LEQ X 0.0) then (ERROR "X must be positive." X) elseif (FGREATERP X 33.0) then (ERROR "X greater than 33." X)) (LET ((COEF (CONSTANT (CL:MAKE-ARRAY 9 (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT) (QUOTE :INITIAL-CONTENTS) (QUOTE (0.03586834 -0.1935278 0.4821994 -0.7567041 0.918207 -0.897057 0.988206 -0.5771917 1.0))))) (FX (FDIFFERENCE X (FIX X))) GLFX) (DECLARE (TYPE FLOATP FX GLFX)) (SETQ GLFX (if (UFGREATERP FX 0.0) then (HORNER FX COEF 8) else 1.0)) (if (FLESSP X 1.0) then (FQUOTIENT GLFX FX) elseif (LEQ X 2.0) then GLFX else (LET ((PROD 1.0) (TERM (FPLUS 1.0 FX)) (TEST (FDIFFERENCE X 0.5))) (DECLARE (TYPE FLOATP PROD TERM TEST)) (do (SETQ PROD (FTIMES PROD TERM)) (SETQ TERM (FPLUS TERM 1.0)) repeatwhile (UFLESSP TERM TEST) finally (RETURN (FTIMES PROD GLFX)))))))
)
(GAMMA-CDF
(LAMBDA (X ETA) (* jop%: " 1-Jun-86 19:13") (* *) (SETQ X (FLOAT X)) (SETQ ETA (FLOAT (OR ETA 1.0))) (if (LEQ X 0.0) then (ERROR "X must be positive." X) elseif (LEQ ETA 0.0) then (ERROR "Eta must be positive." ETA)) (LET* ((FX (FLOAT X)) (FETA ETA) (CETA ETA) (TERM 1.0) (TOTAL 1.0) TEMP) (DECLARE (TYPE FLOATP FX FETA CETA TERM TOTAL TEMP)) (if (AND (UFGEQ FX FETA) (UFGEQ FX 7.0)) then (for I from 1 to (UFIX (if (UFLESSP FX 11.0) then (FPLUS FX FETA -1.0) else (FPLUS FETA 10.0))) do (SETQ CETA (FDIFFERENCE CETA 1.0)) (SETQ TERM (FQUOTIENT (FTIMES TERM CETA) FX)) (SETQ TOTAL (FPLUS TOTAL TERM)) finally (SETQ TERM (FTIMES TERM (FDIFFERENCE CETA 1.0))) (SETQ TOTAL (FPLUS TOTAL (FQUOTIENT TERM (FPLUS FX (UFMINUS CETA) 2.0)))) (RETURN (FDIFFERENCE 1.0 (ANTILOG (SETQ TEMP (FPLUS (UFMINUS FX) (FTIMES (FDIFFERENCE FETA 1.0) (LOG FX)) (LOG TOTAL) (UFMINUS (GAMLOG FETA)))))))) else (repeatwhile (OR (UFGREATERP TEMP 0.5) (UFGREATERP (FQUOTIENT TERM TOTAL) (CONSTANT \IEEE.SHORT.MACHEPS))) do (SETQ CETA (FPLUS CETA 1.0)) (SETQ TEMP (FQUOTIENT FX CETA)) (SETQ TERM (FTIMES TERM TEMP)) (SETQ TOTAL (FPLUS TOTAL TERM)) finally (RETURN (FQUOTIENT (ANTILOG (SETQ TEMP (FPLUS (FTIMES FETA (LOG FX)) (UFMINUS FX) (LOG TOTAL) (UFMINUS (GAMLOG FETA))))) FETA))))))
)
(GAUSSIAN-CDF
(LAMBDA (X MU SIGMA) (* jop%: " 2-Jun-86 22:02") (* *) (SETQ X (FLOAT X)) (SETQ MU (FLOAT (OR MU 0.0))) (SETQ SIGMA (FLOAT (OR SIGMA 1.0))) (if (LEQ SIGMA 0.0) then (ERROR "Sigma must be positive." SIGMA) else (STD-GAUSSIAN-CDF (FQUOTIENT (FDIFFERENCE X MU) SIGMA))))
)
(INV-BETA-CDF
(LAMBDA (PROB A B) (* jop%: " 2-Jun-86 20:43") (* * Algorithm AS 109, Applied Statist. (1977) %, Vol. 26, No. 1) (SETQ PROB (FLOAT PROB)) (SETQ A (FLOAT (OR A 1.0))) (SETQ B (FLOAT (OR B 1.0))) (if (OR (LEQ PROB 0.0) (GEQ PROB 1.0)) then (ERROR "Probability must lie between 0 and 1." PROB) elseif (LEQ A 0.0) then (ERROR "A must be positive." A) elseif (LEQ B 0.0) then (ERROR "B must be positive." B)) (if (OR (FEQP PROB 0.0) (FEQP PROB 1.0) (AND (FEQP A 1.0) (FEQP B 1.0))) then PROB else (LET ((ACU (CONSTANT \IEEE.SHORT.MACHEPS)) (SWAP-P (FGREATERP PROB 0.5)) (BETA (FDIFFERENCE (FPLUS (GAMLOG A) (GAMLOG B)) (GAMLOG (FPLUS A B)))) VALUE) (DECLARE (TYPE FLOATP ACU BETA VALUE)) (* * change tail if necessary) (if SWAP-P then (LET ((AA A)) (SETQ PROB (FDIFFERENCE 1.0 PROB)) (SETQ A B) (SETQ B AA))) (* * Calculate the initial approximation) (LET ((R (SQRT (FMINUS (LOG (FTIMES PROB PROB))))) Y T1 T2 S H W) (DECLARE (TYPE FLOATP R Y T1 T2 S W)) (SETQ Y (FDIFFERENCE R (FQUOTIENT (FPLUS 2.30753 (FTIMES 0.27061 R)) (FPLUS 1.0 (FTIMES (FPLUS 0.99229 (FTIMES 0.04481 R)) R))))) (if (OR (LEQ A 1.0) (LEQ B 1.0)) then (SETQ R (FPLUS B A)) (SETQ T1 (FQUOTIENT 1.0 (FTIMES 9.0 B))) (SETQ T2 (FPLUS 1.0 (UFMINUS T1) (FTIMES Y (SQRT T1)))) (SETQ T1 (FTIMES R T2 T2 T2)) (if (UFLESSP T1 0.0) then (SETQ VALUE (FDIFFERENCE 1.0 (ANTILOG (FQUOTIENT (LOG (FPLUS (FTIMES (FDIFFERENCE 1.0 PROB) B) BETA)) B)))) else (SETQ T1 (FQUOTIENT (FPLUS (TIMES 4.0 A) R -2.0) T1)) (if (UFLESSP T1 1.0) then (SETQ VALUE (ANTILOG (FQUOTIENT (LOG (FPLUS (FTIMES PROB A) BETA)) A))) else (SETQ VALUE (FPLUS 1.0 (UFMINUS (FQUOTIENT 2.0 (FPLUS T1 1.0))))))) else (SETQ R (FQUOTIENT (FPLUS (FTIMES Y Y) -3.0) 6.0)) (SETQ S (FQUOTIENT 1.0 (FPLUS A A -1.0))) (SETQ T1 (FQUOTIENT 1.0 (FPLUS B B -1.0))) (SETQ H (FQUOTIENT 2.0 (FPLUS S T1))) (SETQ W (FDIFFERENCE (FQUOTIENT (FTIMES Y (SQRT (FPLUS H R))) H) (FTIMES (FDIFFERENCE T1 S) (FPLUS R (FQUOTIENT 5.0 6.0) (FQUOTIENT -2.0 (FTIMES 3.0 H)))))) (SETQ VALUE (FQUOTIENT A (FPLUS A (FTIMES B (ANTILOG (FPLUS W W)))))))) (* * Solve for X by a modified newton-raphson method) (if (FLESSP VALUE 1.0E-4) then (SETQ VALUE 1.0E-4) elseif (FGREATERP VALUE 0.9999) then (SETQ VALUE 0.9999)) (bind (R ← (FDIFFERENCE 1.0 A)) (YPREV ← 0.0) (SQ ← 1.0) (PREV ← 1.0) Y OLDVALUE TX declare (TYPE FLOATP R YPREV SQ PREV Y OLDVALUE TX) repeatwhile (AND (UFGEQ PREV ACU) (UFGEQ (FTIMES Y Y) ACU) (NOT (UFEQP TX OLDVALUE))) do (SETQ OLDVALUE VALUE) (SETQ Y (FTIMES (FDIFFERENCE (BETA-CDF VALUE A B) PROB) (ANTILOG (FPLUS BETA (FTIMES R (LOG VALUE)) (FTIMES (FDIFFERENCE 1.0 B) (LOG (FDIFFERENCE 1.0 VALUE))))))) (if (UFLEQ (FTIMES Y YPREV) 0.0) then (SETQ PREV SQ)) (bind (G ← 3.0) ADJ declare (TYPE FLOATP G ADJ) repeatwhile (AND (UFGEQ PREV ACU) (UFGEQ (FTIMES Y Y) PREV) (OR (UFEQP TX 0.0) (UFEQP TX 1.0))) do (repeatwhile (OR (UFGEQ SQ PREV) (UFLESSP TX 0.0) (UFGREATERP TX 1.0)) do (SETQ G (FQUOTIENT G 3.0)) (SETQ ADJ (FTIMES G Y)) (SETQ SQ (FTIMES ADJ ADJ)) (SETQ TX (FDIFFERENCE VALUE ADJ)))) (SETQ VALUE TX) (SETQ YPREV Y)) (if SWAP-P then (FDIFFERENCE 1.0 VALUE) else VALUE))))
)
(INV-CAUCHY-CDF
(LAMBDA (PROB MU SIGMA) (* jop%: "30-May-86 13:21") (* *) (SETQ PROB (FLOAT PROB)) (SETQ MU (FLOAT (OR MU 0.0))) (SETQ SIGMA (FLOAT (OR SIGMA 1.0))) (if (OR (LEQ PROB 0.0) (GEQ PROB 1.0)) then (ERROR "Probability must lie between 0 and 1." PROB) elseif (LEQ SIGMA 0.0) then (ERROR "Sigma must be positive" SIGMA) else (FPLUS MU (FTIMES SIGMA (TAN (FTIMES \IEEE.SHORT.PI (FPLUS PROB -0.5)) T)))))
)
(INV-CHISQUARED-CDF
(LAMBDA (PROB DF) (* jop%: " 1-Jun-86 19:19") (* *) (SETQ PROB (FLOAT PROB)) (SETQ DF (FLOAT (OR DF 1.0))) (if (OR (FLESSP PROB 0.0) (GEQ PROB 1.0)) then (ERROR "Probability must lie between 0 and 1." PROB) elseif (LEQ DF 0.0) then (ERROR "DegreesOfFreedom must be positive." DF) elseif (FEQP PROB 0.0) then 0.0 else (FTIMES 2.0 (INV-GAMMA-CDF PROB (FQUOTIENT DF 2.0)))))
)
(INV-EXPONENTIAL-CDF
(LAMBDA (PROB MU) (* jop%: "25-May-86 16:48") (* *) (SETQ PROB (FLOAT PROB)) (SETQ MU (FLOAT (OR MU 1.0))) (if (OR (FLESSP PROB 0.0) (GEQ PROB 1.0)) then (ERROR "Probability must lie between 0 and 1." PROB) elseif (LEQ MU 0.0) then (ERROR "Mu must be positive." MU) else (FMINUS (FTIMES MU (LOG (FDIFFERENCE 1.0 PROB))))))
)
(INV-F-CDF
(LAMBDA (PROB P Q) (* jop%: " 2-Jun-86 20:42") (* *) (SETQ PROB (FLOAT PROB)) (SETQ P (FLOAT (OR P 1.0))) (SETQ Q (FLOAT (OR Q 1.0))) (if (OR (FLESSP PROB 0.0) (FGREATERP PROB 1.0)) then (ERROR "Probability must lie between 0 and 1." PROB) elseif (LEQ P 0.0) then (ERROR "NumeratorDF must be positive." P) elseif (LEQ Q 0.0) then (ERROR "DenominatorDF must be positive." Q) elseif (FEQP PROB 0.0) then 0.0 elseif (LET ((TEMP (INV-BETA-CDF PROB (FQUOTIENT P 2.0) (FQUOTIENT Q 2.0)))) (FTIMES (FQUOTIENT Q P) (FQUOTIENT TEMP (FDIFFERENCE 1.0 TEMP))))))
)
(INV-GAMMA-CDF
(LAMBDA (PROB ETA) (* jop%: " 2-Jun-86 22:02") (* *) (SETQ PROB (FLOAT PROB)) (SETQ ETA (FLOAT (OR ETA 1.0))) (if (OR (FLESSP PROB 0.0) (GEQ PROB 1.0)) then (ERROR "PROB must lie between 0 and 1") elseif (LEQ ETA 0.0) then (ERROR "Eta must be positive.")) (if (FEQP PROB 0.0) then 0.0 elseif (FEQP ETA 0.5) then (LET ((A (INV-STD-GAUSSIAN-CDF (FTIMES 0.5 (FDIFFERENCE 1.0 PROB))))) (FTIMES A A 0.5)) elseif (FEQP ETA 1.0) then (FMINUS (LOG (FDIFFERENCE 1.0 PROB))) else (LET ((TEMPRESULT (INV-GAMMA-CDF1 PROB ETA))) (if (OR (GEQ ETA 30.0) (AND (GEQ ETA 15.0) (FGREATERP PROB 0.49) (FLESSP PROB 0.99)) (AND (FGREATERP ETA (FDIFFERENCE 25.2 (FTIMES PROB 20.8))) (GEQ PROB 0.01) (LEQ PROB 0.99))) then TEMPRESULT else (INV-GAMMA-CDF2 PROB ETA TEMPRESULT)))))
)
(INV-GAMMA-CDF1
(LAMBDA (PROB ETA) (* jop%: " 2-Jun-86 22:02") (* *) (if (FLESSP ETA 0.5) then (LET* ((A (FQUOTIENT 1.0 (FTIMES 9.0 ETA))) (B (FPLUS 1.0 (FMINUS A) (FTIMES (INV-STD-GAUSSIAN-CDF PROB) (SQRT A)))) (C (FTIMES B B B ETA 2.0))) (if (LEQ C 0.0) then (EXPT PROB (FTIMES 9.0 A)) else C)) else (LET* ((A (FQUOTIENT 0.5 ETA)) (B (FTIMES (SQRT A) (INV-STD-GAUSSIAN-CDF PROB))) (COEFS1 (CONSTANT (CL:MAKE-ARRAY (QUOTE 2) (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT) (QUOTE :INITIAL-CONTENTS) (QUOTE (-0.01425296 0.01264616))))) (COEFS2 (CONSTANT (CL:MAKE-ARRAY (QUOTE 4) (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT) (QUOTE :INITIAL-CONTENTS) (QUOTE (-0.00588609 0.001400483 -0.01091214 -0.002304527))))) (COEFS3 (CONSTANT (CL:MAKE-ARRAY (QUOTE 6) (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT) (QUOTE :INITIAL-CONTENTS) (QUOTE (-2.728484E-4 0.003135411 -0.009699682 0.01316872 0.002618914 -0.2222222))))) (COEFS4 (CONSTANT (CL:MAKE-ARRAY (QUOTE 8) (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT) (QUOTE :INITIAL-CONTENTS) (QUOTE (5.406674E-5 3.483789E-5 -7.274761E-4 0.003292181 -0.008729714 0.0 0.4714045 1.0))))) (C (FPLUS (FTIMES (FPLUS (FTIMES (FPLUS (FTIMES (HORNER B COEFS1 1) A) (HORNER B COEFS2 3)) A) (HORNER B COEFS3 5)) A) (HORNER B COEFS4 7)))) (FTIMES C C C ETA))))
)
(INV-GAMMA-CDF2
(LAMBDA (PROB ETA TEMPRESULT) (* jop%: "23-Jun-86 13:30") (* *) (LET* ((GAMMATEMP (GAMMA-CDF TEMPRESULT ETA)) (PDIFF (FDIFFERENCE PROB GAMMATEMP)) (ABSPDIFF (ABS PDIFF)) (OLDABSCORRECTION MIN.FLOAT) (OLDTEMP TEMPRESULT) CORRECTION ABSCORRECTION) (for I from 1 to 50 repeatwhile (AND (FGREATERP ABSPDIFF 1.0E-8) (FGREATERP OLDABSCORRECTION (FTIMES OLDTEMP 1.0E-6))) do (if (AND (GEQ ABSPDIFF GAMMATEMP) (NOT (EQ I 1))) then (SETQ CORRECTION (FMINUS (FQUOTIENT CORRECTION 2.0))) else (SETQ CORRECTION (FPLUS (LOG ABSPDIFF) (FMINUS (FTIMES (FDIFFERENCE ETA 1.0) (LOG TEMPRESULT))) TEMPRESULT (GAMLOG ETA))) (SETQ CORRECTION (if (FLESSP PDIFF 0.0) then (FMINUS (ANTILOG CORRECTION)) else (ANTILOG CORRECTION))) (if (NOT (EQ (FLESSP PDIFF 0.0) (FLESSP CORRECTION 0.0))) then (SETQ CORRECTION (if (FLESSP PDIFF 0.0) then (FMINUS (FQUOTIENT ABSCORRECTION 2.0)) else (FQUOTIENT ABSCORRECTION 2.0))))) (SETQ ABSCORRECTION (ABS CORRECTION)) (SETQ OLDABSCORRECTION ABSCORRECTION) (SETQ GAMMATEMP ABSPDIFF) (SETQ TEMPRESULT (FMAX (FPLUS TEMPRESULT CORRECTION) \IEEE.SHORT.LEAST.NORMAL)) (SETQ OLDTEMP TEMPRESULT) (SETQ GAMMATEMP (GAMMA-CDF TEMPRESULT ETA)) (SETQ PDIFF (FDIFFERENCE PROB GAMMATEMP)) (SETQ ABSPDIFF (ABS PDIFF)) finally (RETURN OLDTEMP))))
)
(INV-GAUSSIAN-CDF
(LAMBDA (PROB MU SIGMA) (* jop%: " 2-Jun-86 22:02") (* *) (SETQ PROB (FLOAT PROB)) (SETQ MU (FLOAT (OR MU 0.0))) (SETQ SIGMA (FLOAT (OR SIGMA 1.0))) (if (OR (LEQ PROB 0.0) (GEQ PROB 1.0)) then (ERROR "Probability must lie between 0 and 1." PROB) elseif (LEQ SIGMA 0.0) then (ERROR "Scale must be positive." SIGMA) else (FPLUS MU (FTIMES SIGMA (INV-STD-GAUSSIAN-CDF PROB)))))
)
(INV-LOGISTIC-CDF
(LAMBDA (PROB MU SIGMA) (* jop%: "30-May-86 13:24") (* *) (SETQ PROB (FLOAT PROB)) (SETQ MU (FLOAT (OR MU 0.0))) (SETQ SIGMA (FLOAT (OR SIGMA 1.0))) (if (OR (LEQ PROB 0.0) (GEQ PROB 1.0)) then (ERROR "Probability must lie between 0 and 1." PROB) elseif (LEQ SIGMA 0.0) then (ERROR "Sigma must be positive." SIGMA) else (FPLUS MU (FTIMES SIGMA (LOG (FQUOTIENT PROB (FDIFFERENCE 1.0 PROB)))))))
)
(INV-LOGNORMAL-CDF
(LAMBDA (PROB MU SIGMA) (* jop%: " 2-Jun-86 22:02") (* *) (SETQ PROB (FLOAT PROB)) (SETQ MU (FLOAT (OR MU 0.0))) (SETQ SIGMA (FLOAT (OR SIGMA 1.0))) (if (OR (LEQ PROB 0.0) (GEQ PROB 1.0)) then (ERROR "Probability must lie between 0 and 1." PROB) elseif (LEQ SIGMA 0.0) then (ERROR "Sigma must be positive." SIGMA) else (ANTILOG (FPLUS MU (FTIMES SIGMA (INV-STD-GAUSSIAN-CDF PROB))))))
)
(INV-STD-GAUSSIAN-CDF
(LAMBDA (P) (* jop%: " 3-Jun-86 15:13") (* * Algorithm AS 11, Applied Stat. (1977) %, Vol. 26, N0. 1) (* * produces standard normal deviate corresponding to lower tail area of p.) (SETQ P (FLOAT P)) (LET ((N1 (CONSTANT (CL:MAKE-ARRAY 4 (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT) (QUOTE :INITIAL-CONTENTS) (QUOTE (-25.44106 41.3912 -18.615 2.506628))))) (D1 (CONSTANT (CL:MAKE-ARRAY 5 (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT) (QUOTE :INITIAL-CONTENTS) (QUOTE (3.130829 -21.06224 23.08337 -8.473512 1.0))))) (N2 (CONSTANT (CL:MAKE-ARRAY 4 (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT) (QUOTE :INITIAL-CONTENTS) (QUOTE (2.321213 4.850141 -2.297965 -2.78719))))) (D2 (CONSTANT (CL:MAKE-ARRAY 3 (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT) (QUOTE :INITIAL-CONTENTS) (QUOTE (1.637068 3.543889 1.0))))) (Q (FDIFFERENCE P 0.5)) R VALUE) (DECLARE (TYPE FLOATP Q R VALUE)) (if (UFLESSP (UFABS Q) 0.42) then (SETQ R (FTIMES Q Q)) (SETQ VALUE (FQUOTIENT (FTIMES Q (HORNER R N1 3)) (HORNER R D1 4))) else (SETQ R (SQRT (FMINUS (LOG (if (UFGREATERP Q 0.0) then (FDIFFERENCE 1.0 P) else P))))) (SETQ VALUE (FQUOTIENT (HORNER R N2 3) (HORNER R D2 2))) (if (UFLESSP Q 0.0) then (UFMINUS VALUE) else VALUE))))
)
(INV-T-CDF
(LAMBDA (PROB DF) (* jop%: " 6-Jun-86 16:07") (* * Adapted from the S routine qtdis) (SETQ PROB (FLOAT PROB)) (SETQ DF (FLOAT (OR DF 1.0))) (if (OR (LEQ PROB 0.0) (GEQ PROB 1.0)) then (ERROR "Probability must lie between 0 and 1." PROB) elseif (FLESSP DF 1.0) then (ERROR "DF must be greater than 1.0" DF)) (LET ((HALFPI 1.570796) (ZROTOL \IEEE.SHORT.LEAST.NORMAL) (FN DF) (PR (FTIMES 2.0 (FDIFFERENCE 1.0 PROB))) VALUE) (DECLARE (TYPE FLOATP HALFPI ZROTOL FN PR VALUE)) (if (FLESSP PROB 0.5) then (SETQ PR (FTIMES PROB 2.0))) (if (UFLEQ (UFABS (FDIFFERENCE FN 2.0)) ZROTOL) then (* For DF = 2) (SETQ VALUE (SQRT (SETQ VALUE (FDIFFERENCE (FQUOTIENT 2.0 (FTIMES PR (FDIFFERENCE 2.0 PR))) 2.0)))) elseif (UFLEQ (UFABS (FDIFFERENCE FN 1.0)) ZROTOL) then (* For DF = 1) (SETQ PR (FTIMES PR HALFPI)) (SETQ VALUE (FQUOTIENT (COS PR T) (SIN PR T))) else (LET ((A (FQUOTIENT 1.0 (FDIFFERENCE FN 0.5))) B C D X Y PR2 FNP2) (DECLARE (TYPE FLOATP A B C D X Y PR2 FNP2)) (SETQ B (FQUOTIENT 48.0 (FTIMES A A))) (SETQ C (FPLUS (FTIMES (FPLUS (FTIMES (FPLUS (FQUOTIENT (FTIMES 20700.0 A) B) -98.0) A) -16.0) A) 96.36)) (SETQ D (FTIMES (FPLUS (FQUOTIENT (FPLUS (FQUOTIENT 94.5 (FPLUS B C)) -3.0) B) 1.0) (SQRT (FTIMES A HALFPI)) FN)) (SETQ X (FTIMES D PR)) (SETQ Y (EXPT X (FQUOTIENT 2.0 FN))) (if (UFLEQ Y (FPLUS A 0.05)) then (SETQ FNP2 (FPLUS FN 2.0)) (SETQ Y (FPLUS (FQUOTIENT (FTIMES (FPLUS (FTIMES (FPLUS (FQUOTIENT 1.0 (FTIMES (FPLUS (FQUOTIENT (FPLUS FN 6.0) (FPLUS FN Y)) (FTIMES -0.089 D) -0.8220001) FNP2 3.0)) (FQUOTIENT 0.5 (FPLUS FN 4.0))) Y) -1.0) (FPLUS FN 1.0)) FNP2) (FQUOTIENT 1.0 Y))) else (* asymptotic inverse expansion about normal) (SETQ PR2 (FTIMES 0.5 PR)) (SETQ X (INV-STD-GAUSSIAN-CDF PR2)) (SETQ Y (FTIMES X X)) (if (UFLESSP FN 5.0) then (SETQ C (FPLUS C (FTIMES 0.3 (FDIFFERENCE FN 4.5) (FPLUS X 0.6))))) (SETQ C (FPLUS (FTIMES (FPLUS (FTIMES (FPLUS (FTIMES (FPLUS (FTIMES 0.05 D X) -5.0) X) -7.0) X) -2.0) X) B C)) (SETQ Y (FTIMES (FPLUS (FQUOTIENT (FPLUS (FQUOTIENT (FPLUS (FTIMES (FPLUS (FTIMES (FPLUS (FTIMES 0.4 Y) 6.3) Y) 36.0) Y) 94.5) C) (UFMINUS Y) -3.0) B) 1.0) X)) (SETQ Y (FTIMES A Y Y)) (if (UFLEQ Y 0.002) then (SETQ Y (FPLUS Y (FTIMES 0.5 Y Y))) else (SETQ Y (FDIFFERENCE (ANTILOG Y) 1.0)))) (SETQ VALUE (SQRT (FTIMES FN Y))))) (if (FLESSP PROB 0.5) then (UFMINUS VALUE) else VALUE)))
)
(INV-UNIFORM-CDF
(LAMBDA (PROB LOW HIGH) (* jop%: "25-May-86 15:43") (* *) (SETQ PROB (FLOAT PROB)) (SETQ LOW (FLOAT (OR LOW 0.0))) (SETQ HIGH (FLOAT (OR HIGH 1.0))) (if (OR (FLESSP PROB 0.0) (FGREATERP PROB 1.0)) then (ERROR "Probability must lie between 0 and 1.") elseif (LEQ HIGH LOW) then (ERROR "No effective calculation range.") else (FPLUS (FTIMES PROB (FDIFFERENCE HIGH LOW)) LOW)))
)
(LOGISTIC-CDF
(LAMBDA (X MU SIGMA) (* jop%: "30-May-86 13:07") (* *) (SETQ X (FLOAT X)) (SETQ MU (FLOAT (OR MU 0.0))) (SETQ SIGMA (FLOAT (OR SIGMA 1.0))) (if (LEQ SIGMA 0.0) then (ERROR "Sigma must be positive." SIGMA) else (FQUOTIENT 1.0 (FPLUS 1.0 (ANTILOG (FQUOTIENT (FDIFFERENCE MU X) SIGMA))))))
)
(LOGNORMAL-CDF
(LAMBDA (X MU SIGMA) (* jop%: " 2-Jun-86 22:02") (* *) (SETQ X (FLOAT X)) (SETQ MU (FLOAT (OR MU 0.0))) (SETQ SIGMA (FLOAT (OR SIGMA 1.0))) (if (FLESSP X 0.0) then (ERROR "X must be positive." X) elseif (LEQ SIGMA 0.0) then (ERROR "Sigma must be positive." SIGMA) elseif (FEQP X 0.0) then 0.0 else (STD-GAUSSIAN-CDF (FQUOTIENT (FDIFFERENCE (LOG X) MU) SIGMA))))
)
(STD-GAUSSIAN-CDF
(LAMBDA (X) (* jop%: " 3-Jun-86 15:05") (* * Algorithm taken from Appl. Statistics (1973) Vol. 22, No.3) (* * Evaluates the tail area of the standard gaussian density from -infinity to x) (SETQ X (FLOAT X)) (* 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 (FABS X)) (CHANGESIGN (FLESSP X 0.0)) FY VALUE) (DECLARE (TYPE FLOATP FX FY VALUE)) (SETQ VALUE (if (OR (UFGREATERP FX UTZERO) (AND CHANGESIGN (UFGREATERP FX LTONE))) then 0.0 else (SETQ FY (FTIMES 0.5 FX FX)) (if (UFLESSP FX 1.28) then (FPLUS 0.5 (FTIMES (UFMINUS FX) (FPLUS 0.3989423 (FQUOTIENT (FTIMES -0.3999035 FY) (FPLUS FY 5.758855 (FQUOTIENT -29.82136 (FPLUS FY 2.624331 (FQUOTIENT 48.696 (FPLUS FY 5.928857))))))))) else (FQUOTIENT (FTIMES 0.3989423 (ANTILOG (FMINUS FY))) (FPLUS FX -3.8052E-8 (FQUOTIENT 1.000006 (FPLUS FX 3.980648E-4 (FQUOTIENT 1.986154 (FPLUS FX -0.1516791 (FQUOTIENT 5.293304 (FPLUS FX 4.838592 (FQUOTIENT -15.1509 (FPLUS FX 0.7423809 (FQUOTIENT 30.78993 (FPLUS FX 3.990194))))))))))))))) (if CHANGESIGN then VALUE else (FDIFFERENCE 1.0 VALUE))))
)
(T-CDF
(LAMBDA (X DF) (* jop%: "23-Jun-86 13:31") (* * adaptation of ACM algorithm 395 (G. W. Hill, v. 13 (1970) %, |617-619|)) (SETQ X (FLOAT X)) (SETQ DF (FLOAT (OR DF 1.0))) (if (LEQ DF 0.0) then (ERROR "DFmust be greater than zero" DF)) (LET ((ZROTOL (CONSTANT \IEEE.SHORT.LEAST.NORMAL)) (A1PLUS (CONSTANT (FPLUS 1.0 \IEEE.SHORT.MACHEPS))) (TWO-OVER-PI (CONSTANT (FQUOTIENT 2.0 \IEEE.SHORT.PI))) (TP X) (TVALUE X) (N (FIX DF)) (FN DF) (Z 1.0) Y B A STUDEN) (DECLARE (TYPE FLOATP ZROTOL A1PLUS TWO-OVER-PI TP TVALUE FN Z Y B A STUDEN)) (SETQ TVALUE (FTIMES TVALUE TVALUE)) (SETQ Y (FQUOTIENT TVALUE FN)) (SETQ B (FPLUS 1.0 Y)) (if (OR (UFGEQ (UFABS (FDIFFERENCE FN (FLOAT N))) 1.0E-6) (AND (IGEQ N 20) (UFLESSP TVALUE FN)) (IGREATERP N 200)) then (* asymptotic series for large or non-integer n) (if (UFGREATERP Y 1.0E-6) then (SETQ Y (LOG B))) (SETQ A (FDIFFERENCE FN 0.5)) (SETQ B (FTIMES 48.0 A A)) (SETQ Y (FTIMES A Y)) (SETQ Y (FTIMES (FPLUS (FQUOTIENT (FPLUS (FQUOTIENT (FPLUS (FTIMES (FPLUS (FTIMES (FPLUS (FTIMES -0.4 Y) -3.3) Y) -24.0) Y) -85.5) (FPLUS (FTIMES 0.8 Y Y) 100.0 B)) Y 3.0) B) 1.0) (SQRT Y))) (SETQ STUDEN (FTIMES 2.0 (STD-GAUSSIAN-CDF (FMINUS Y)))) else (if (AND (ILESSP N 20) (UFLESSP TVALUE 4.0)) then (* nested summation of "cosine" series) (SETQ A (SQRT Y)) (SETQ Y A) (if (EQ N 1) then (SETQ A 0.0)) else (* "tail" series expansion for large t-values) (SETQ A (SQRT B)) (SETQ Y (FTIMES A FN)) (for J from 2 by 2 until (UFLEQ (UFABS (FDIFFERENCE A Z)) ZROTOL) do (SETQ Z A) (SETQ Y (FQUOTIENT (FTIMES Y (FLOAT (SUB1 J))) (FTIMES B (FLOAT J)))) (SETQ A (FPLUS A (FQUOTIENT Y (FPLUS FN (FLOAT J)))))) (SETQ FN (FPLUS FN 2.0)) (SETQ Y 0.0) (SETQ Z 0.0) (SETQ A (UFMINUS A))) (until (UFLEQ (SETQ FN (FDIFFERENCE FN 2.0)) A1PLUS) do (SETQ A (FPLUS Y (FTIMES A (FQUOTIENT (FDIFFERENCE FN 1.0) (FTIMES B FN)))))) (if (UFGREATERP (FABS FN) ZROTOL) then (SETQ A (FTIMES TWO-OVER-PI (FPLUS (ARCTAN Y T) (FQUOTIENT A B)))) else (SETQ A (FQUOTIENT A (SQRT B)))) (SETQ STUDEN (FDIFFERENCE Z A))) (if (UFGEQ TP 0.0) then (FDIFFERENCE 1.0 (FTIMES 0.5 STUDEN)) else (FTIMES 0.5 STUDEN))))
)
(UNIFORM-CDF
(LAMBDA (X LOW HIGH) (* jop%: " 2-Jun-86 20:45") (* *) (SETQ X (FLOAT X)) (SETQ LOW (FLOAT (OR LOW 0.0))) (SETQ HIGH (FLOAT (OR HIGH 1.0))) (if (LEQ HIGH LOW) then (ERROR "No effective calculation range.") elseif (OR (LEQ X LOW) (GEQ X HIGH)) then 0.0 else (FQUOTIENT (FDIFFERENCE X LOW) (FDIFFERENCE HIGH LOW))))
)
)
(DECLARE%: EVAL@COMPILE
(PUTPROPS HORNER DMACRO ((X POLY CL::DEGREE) (\FLOATBOX ((OPCODES UBFLOAT3 0) (\FLOATUNBOX X) (%%ARRAY-BASE POLY) CL::DEGREE))))
)
(DECLARE%: EVAL@COMPILE
(RPAQ \IEEE.SHORT.MACHEPS (EXPT 2.0 -23))
(RPAQ \IEEE.SHORT.LEAST.NORMAL (EXPT 2.0 -126))
(RPAQQ \IEEE.SHORT.BASE 2)
(RPAQQ \IEEE.SHORT.MANTISSA 23)
(RPAQQ \IEEE.SHORT.EXPONENT 8)
(RPAQQ \IEEE.SHORT.LEAST.EXPONENT -126)
(RPAQQ \IEEE.SHORT.GREATEST.EXPONENT 127)
(RPAQQ \IEEE.SHORT.EXPONENT.BIAS 127)
(RPAQQ \IEEE.SHORT.PI 3.141593)
(CONSTANTS (\IEEE.SHORT.MACHEPS (EXPT 2.0 -23)) (\IEEE.SHORT.LEAST.NORMAL (EXPT 2.0 -126)) (\IEEE.SHORT.BASE 2) (\IEEE.SHORT.MANTISSA 23) (\IEEE.SHORT.EXPONENT 8) (\IEEE.SHORT.LEAST.EXPONENT -126) (\IEEE.SHORT.GREATEST.EXPONENT 127) (\IEEE.SHORT.EXPONENT.BIAS 127) (\IEEE.SHORT.PI 3.141593))
)
(DECLARE%: DONTEVAL@LOAD DOEVAL@COMPILE DONTCOPY
(DECLARE%: DOEVAL@COMPILE DONTCOPY
(LOCALVARS . T)
)
)
(PUTPROPS CDFS COPYRIGHT ("Xerox Corporation" 1986 1988))
(DECLARE%: DONTCOPY
(FILEMAP (NIL (1367 25813 (BETA-CDF 1377 . 3871) (CAUCHY-CDF 3873 . 4178) (CHISQUARED-CDF 4180 . 4511)
(EXPONENTIAL-CDF 4513 . 4808) (F-CDF 4810 . 5292) (GAMLOG 5294 . 5723) (GAMMA 5725 . 6576) (GAMMA-CDF
6578 . 7853) (GAUSSIAN-CDF 7855 . 8142) (INV-BETA-CDF 8144 . 11250) (INV-CAUCHY-CDF 11252 . 11669) (
INV-CHISQUARED-CDF 11671 . 12068) (INV-EXPONENTIAL-CDF 12070 . 12419) (INV-F-CDF 12421 . 12988) (
INV-GAMMA-CDF 12990 . 13766) (INV-GAMMA-CDF1 13768 . 15016) (INV-GAMMA-CDF2 15018 . 16282) (
INV-GAUSSIAN-CDF 16284 . 16682) (INV-LOGISTIC-CDF 16684 . 17100) (INV-LOGNORMAL-CDF 17102 . 17511) (
INV-STD-GAUSSIAN-CDF 17513 . 18701) (INV-T-CDF 18703 . 21033) (INV-UNIFORM-CDF 21035 . 21432) (
LOGISTIC-CDF 21434 . 21740) (LOGNORMAL-CDF 21742 . 22124) (STD-GAUSSIAN-CDF 22126 . 23367) (T-CDF
23369 . 25477) (UNIFORM-CDF 25479 . 25811)))))
STOP