----- Original Message -----
From: Camilo La Rota <[EMAIL PROTECTED]>
To: David A. Heiser <[EMAIL PROTECTED]>
Sent: Tuesday, November 23, 1999 7:41 PM
Subject: Re: quantiles, how to calculate?
> Hi David,
>
> Thank you for your answer.
>
> Yes of course I am interested.
>
> Can you also give me a reference or tell me how this is done?
>
> Thank you very much,
>
> Camilo
>
__________________________________________
Attached is a text file listing the two subroutines involved. I do all my
calcs in double precision which explains the # signs.
DAH
DEFINT I-N
DEFLNG H
DEFSTR F, S
DEFDBL A-E, O-R, T-Z
SUB NormDist (PV, KK, ZV, DV, KV)
REM Inverse normal distribution, Input, PV, Output ZV (standarized normal deviate), DV
REM PV, input 0-1
REM KK, control, =0 approximate ZV, no DV
REM =-1, precise value, DV calculated
REM ZV, deviate or sigma value -1E+38 TO 1E+38
REM DV distribution density
REM KV, =-1 if correct input
REM KV=0, error in PV input
SELECT CASE PV
CASE IS < 0#
ZV = -1E+38 '1213'
DV = 0
PRINT "PM is less than ZERO"
KV = 0
EXIT SUB
CASE IS > 1#
ZV = 1E+38
DV = 0
PRINT "PM greater than one"
KV = 0
EXIT SUB
CASE IS = 0#
ZV = -1E+38
DV = 0#
KV = -1
EXIT SUB
CASE IS = .5#
ZV = 0#
DV = .39894228040143#
KV = -1
EXIT SUB
CASE IS = 1#
ZV = 1E+38
DV = 0#
KV = -1
EXIT SUB
CASE IS < .5#
T0 = PV
T9 = -1#
GOTO inverse
CASE IS > .5#
T0 = 1# - PV
T9 = 1#
GOTO inverse
END SELECT
inverse:
REM approximation
T1 = SQR(LOG(1# / (T0 * T0)))
T2 = 2.515517# + T1 * (.802853# + T1 * .010328#)
T3 = 1# + T1 * (1.432788# + T1 * (.189269# + T1 * .001308#))
T4 = T1 - T2 / T3
ZV = T9 * T4
DV = 0#
IF NOT KK THEN
KV = -1
EXIT SUB
ELSE
REM Newton's method to close in on precise ZV value
DO
CALL NormDistPrb(ZV, PVC, QVC, DV)
DP = ABS(PVC - PV)
REM accuracy can be changed by resetting the DP test value
IF DP < .00000000000001# THEN
KV = -1
EXIT DO
ELSE
IF DV = 0# THEN
PRINT "Error in routine"
PRINT "Current Z value"; ZV
STOP
KV = 0
EXIT SUB
ELSE
ZV = ZV - (PVC - PV) / DV
END IF
END IF
LOOP
END IF
END SUB
SUB NormDistPrb (ZD, PV, QV, ZI)
REM Normal distribution, N(0,1), Input ZD, Output ZI, PV and QV
DIM AT(4), BT(4), AC(2), BC(2)
X0 = ABS(ZD)
IF X0 > 17# THEN
ZI = 0#
IF ZD < 0# THEN
PV = 0#
QV = 1#
EXIT SUB
ELSE
PV = 1#
QV = 0#
EXIT SUB
END IF
ELSE
REM probability density, ZI
XA = -ZD * ZD / 2!
IF X0 > 3.6 THEN
ZI = .39894228040142# * EXP(XA)
ELSE
TN = 0#
T1 = 1#
EX = 1#
DO
TN = TN + 1#
T1 = T1 * XA / TN
EX = EX + T1
IF (ABS(T1 / EX)) < 1E-14 THEN EXIT DO
LOOP
ZI = .39894228040142# * EX
END IF
IF X0 < 2# THEN
REM normal series
NN = 0
T1 = ZI * ZD
PV = .5# + T1
DO
NN = NN + 1
T1 = T1 * ZD * ZD / (2# * CDBL(NN) + 1#)
PV = PV + T1
IF (ABS(T1)) < 1E-15 THEN EXIT DO
LOOP
QV = 1# - PV
EXIT SUB
ELSE
AT(2) = 1#: AT(1) = 0#
BT(2) = X0: BT(1) = 1#
BC(1) = X0: BC(2) = X0
NN = 0: AO = 0#
DO
REM continuing fraction
NN = NN + 1
AC(2) = CDBL(2 * NN)
AC(1) = AC(2) - 1#
REM computes continuing fraction sums
AT(3) = BC(1) * AT(2) + AC(1) * AT(1)
BT(3) = BC(1) * BT(2) + AC(1) * BT(1)
AT(4) = BC(2) * AT(3) + AC(2) * AT(2)
BT(4) = BC(2) * BT(3) + AC(2) * BT(2)
AT(1) = AT(3) / BT(4)
BT(1) = BT(3) / BT(4)
AT(2) = AT(4) / BT(4)
BT(2) = 1#
DF = ABS((AT(2) - AO) / AT(2))
AO = AT(2)
IF DF < 1E-14 THEN EXIT DO
LOOP
IF ZD < 0# THEN
PV = ZI * AT(2)
QV = 1# - PV
EXIT SUB
ELSE
QV = ZI * AT(2)
PV = 1# - QV
EXIT SUB
END IF
END IF
END IF
END SUB
