(FILECREATED " 7-Oct-85 22:09:17" {QV}<IDL>SOURCES>DIST.;9 7864
changes to: (FNS RANDN FPROB.LISP TPROB)
previous date: " 5-Sep-85 18:24:58" {QV}<IDL>SOURCES>DIST.;8)
(* Copyright (c) 1983, 1984, 1985 by Xerox Corporation. All rights reserved.)
(PRETTYCOMPRINT DISTCOMS)
(RPAQQ DISTCOMS ((* Distribution and generator functions)
(* Distribution functions rewritten from code by Rulifson from algorithms in CACM)
(FNS FPROB FPROB.LISP NPROB TPROB TVALUE)
(* Generator functions for non uniform distributions)
(FNS RANDN)))
(* Distribution and generator functions)
(* Distribution functions rewritten from code by Rulifson from algorithms in CACM)
(DEFINEQ
(FPROB
[ULAMBDA ((F (EXPECTS SCALAR))
(DFNUM (EXPECTS SCALAR))
(DFDEN (EXPECTS SCALAR))
(RETURNS SCALAR))
(* rmk: "16-AUG-78 00:52")
(AND F DFNUM DFDEN (FGREATERP F 0.0)
(FGREATERP DFNUM 0.0)
(FGREATERP DFDEN 0.0)
(FPROB.LISP F DFNUM DFDEN))])
(FPROB.LISP
[DLAMBDA ((F ARITH)
(DFNUM ARITH)
(DFDEN ARITH))
(* jop: " 5-Sep-85 17:12")
(* CACM algorithm 322, Feb 68 pg 116, remark Jan 69 pg 39, and remark Feb 71 pg 117 -
For Student's T-distribution call (FPROB.LISP T*T 1 DFDEN) -
The .LISP version is explicit for calls from system, e.g. from ANOVA and TPROB)
(DPROG ((A (IPLUS (ITIMES (IQUOTIENT DFNUM 2)
2)
(IMINUS DFNUM)
2) INTEGER)
(B (IPLUS (ITIMES (IQUOTIENT DFDEN 2)
2)
(IMINUS DFDEN)
2) INTEGER)
(I NIL INTEGER)
(J NIL INTEGER)
(W (FTIMES F (FQUOTIENT DFNUM DFDEN)) FLOATING)
(Y NIL FLOATING)
(Z NIL FLOATING)
(D NIL FLOATING)
(P NIL FLOATING)
(ZK NIL FLOATING))
(SETQ Z (FQUOTIENT 1.0 (FPLUS 1.0 W)))
[if (EQ A 1)
then [if (EQ B 1)
then (SETQ P (SQRT W))
(SETQ Y .3183099)
(SETQ D (FTIMES Y (FQUOTIENT Z P)))
(SETQ P (FTIMES 2.0 Y (ARCTAN P T)))
else (SETQ P (SQRT (FTIMES W Z)))
(SETQ D (FTIMES .5 P (FQUOTIENT Z W]
else (if (EQ B 1)
then (SETQ P (SQRT Z))
(SETQ D (FTIMES .5 Z P))
(SETQ P (FDIFFERENCE 1.0 P))
else (SETQ D (FTIMES Z Z))
(SETQ P (FTIMES W Z]
(SETQ Y (FTIMES 2.0 (FQUOTIENT W Z)))
[if (EQ A 1)
then [for J from (IPLUS B 2) by 2 to DFDEN
do (SETQ D (FTIMES (FPLUS 1 (FQUOTIENT A (FDIFFERENCE J 2.0)))
D Z))
(SETQ P (FPLUS P (FTIMES D (FQUOTIENT Y (SUB1 J]
else (SETQ ZK (EXPT Z (IQUOTIENT (SUB1 DFDEN)
2)))
(SETQ D (FTIMES D ZK (FQUOTIENT DFDEN B)))
(SETQ P (FPLUS (FTIMES P ZK)
(FTIMES W Z (FQUOTIENT (FDIFFERENCE ZK 1.0)
(FDIFFERENCE Z 1.0]
(SETQ Y (FTIMES W Z))
(SETQ Z (FQUOTIENT 2.0 Z))
(SETQ B (IDIFFERENCE DFDEN 2))
[for I from (IPLUS A 2) by 2 to DFNUM
do (SETQ J (IPLUS I B))
[SETQ D (FTIMES Y D (FQUOTIENT J (IDIFFERENCE I 2]
(SETQ P (FDIFFERENCE P (FTIMES Z (FQUOTIENT D J]
(RETURN (if (FGTP P 1.0)
then 0.0
elseif (FGTP 0.0 P)
then 1.0
else (FDIFFERENCE 1.0 P))))])
(NPROB
[ULAMBDA ((Z (EXPECTS SCALAR)))
(* bas: "11-FEB-83 11:14")
(* One-tailed normal probability function, Handbook of
Mathematical Functions, Abramowitz and Stegun, 26.2.19,
pg. 932)
(AND
Z
(DPROGN ((Z ARITH))
(FQUOTIENT
(EXPT
(FPLUS
(FTIMES (FPLUS (FTIMES (FPLUS (FTIMES (FPLUS (FTIMES
(FPLUS (FTIMES (FPLUS (FTIMES 5.383E-6
Z)
.0000488906)
Z)
.0000380036)
Z)
.003277626)
Z)
.02114101)
Z)
.04986735)
Z)
1.0)
-16)
2.0)))])
(TPROB
[ULAMBDA ((X (EXPECTS SCALAR))
(DF (EXPECTS SCALAR)))
(* jop: " 5-Sep-85 17:20")
(* Probability from Student's T distribution.)
(AND X DF (GREATERP DF 0)
(FPROB.LISP (FTIMES X X)
1 DF))])
(TVALUE
[ULAMBDA ((P (EXPECTS SCALAR))
(DF (EXPECTS SCALAR) SPECIAL))
(* bas: "11-FEB-83 11:15")
(* Computes X, where (TPROB X DF) = P)
(AND P DF (GREATERP DF 0)
(FNINVERSE (FUNCTION [LAMBDA (X)
(TPROB X DF])
P 2))])
)
(* Generator functions for non uniform distributions)
(DEFINEQ
(RANDN
[ULAMBDA ((MEAN (EXPECTS SCALAR))
(STDEV (EXPECTS SCALAR)))
(* jop: " 5-Sep-85 17:21")
(* Polar method of generating random normal deviates from Knuth Vol 2 p 104 Made a little bit more hairy by reuse of
boxes from SQRT RAND and LOG whenever possible)
(DPROG ((SAVE (CONSTANT (LIST NIL)) LISTP (* Location for saving the second result))
THEN (VAL (OR (PROG1 (CAR SAVE)
(FRPLACA SAVE))
(DPROG ((V1 NIL (SATISFIES (FLOATP V1))
(* Suppressing Ron K))
(V2 (RAND -1.0 1.0) (SATISFIES (FLOATP V2))
(* Ditto))
(S NIL FLOATP)
(RETURNS FLOATP))
(do (SETQ V1 V2) (* V1 and V2 are decled funny to keep SETQ macros from
allocating a static boxes for them)
(SETQ V2 (RAND -1.0 1.0))
(SETQ S (FPLUS (FTIMES V1 V1)
(FTIMES V2 V2)))
repeatuntil (FGTP 1.0 S))
(SETQ S (SQRT (FQUOTIENT (FTIMES -2.0 (LOG S))
S)))
(replace F of V1 with (FTIMES V1 S))
(* These are OK because V1 and V2 are always boxes from
RAND)
(FRPLACA SAVE V1) (* Save the second value away for the next time)
(replace F of V2 with (FTIMES V2 S))
(RETURN V2))) (SATISFIES (FLOATP VAL))
(* VAL is bound to a box from RAND always))
(RETURNS FLOATP))
(replace F of VAL with (FPLUS (OR MEAN 0.0)
(FTIMES (OR STDEV 1.0)
VAL)))
(RETURN VAL))])
)
(PUTPROPS DIST COPYRIGHT ("Xerox Corporation" 1983 1984 1985))
(DECLARE: DONTCOPY
(FILEMAP (NIL (706 5506 (FPROB 716 . 1096) (FPROB.LISP 1098 . 3874) (NPROB 3876 . 4734) (TPROB 4736 .
5107) (TVALUE 5109 . 5504)) (5569 7779 (RANDN 5579 . 7777)))))
STOP