'********************************************************************* '* This program calculates the standardized normal law probabilities * '* (mean=0, standard deviation=1) for -37.0 <= X <= 37.0 * '* ----------------------------------------------------------------- * '* SAMPLE RUN: * '* X=5. Number of standard deviations * '* P1 Probability to the left of X * '* P2 Probability to the right of X * '* Q Probability Density for X * '* (Use successively function alnorm and subroutines NORMP and * '* NPROB). * '* * '* P1= .9999997133484282 * '* P2= 2.866515718673977D-07 * '* * '* P1= .9999997133484281 * '* P2= 2.866515718923335D-07 * '* Q= 1.486719514734297D-06 * '* * '* P1= .9999997133484282 * '* P2= 2.866515718673977D-07 * '* Q= 1.486719514673059D-06 * '* * '* ----------------------------------------------------------------- * '* Ref.: Journal of Applied Statistics (1973) vol22 no.3. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '********************************************************************* 'Program TNormal DEFDBL A-H, O-Z DEFINT I-N zero = 0# one = 1# half = .5 ltone = 7# utzero = 18.66 con = 1.28 ' P1, P2, Q, X: Double ' IUP: Integer (0 or 1) ' This file includes the Applied Statistics algorithm AS 66 for calculating ' the tail area under the normal curve, and two alternative routines which ' give higher accuracy. The latter have been contributed by Alan Miller of ' CSIRO Division of Mathematics & Statistics, Clayton, Victoria. Notice ' that each function or routine has different call arguments. X = 5# IUP = 0 GOSUB 500 'P1= alnorm(X,IUP) P1 = alnorm IUP = 1 GOSUB 500 'P2= alnorm(X,IUP) P2 = alnorm CLS PRINT PRINT " P1="; P1 PRINT " P2="; P2 GOSUB 1000 'call NORMP(X, P, Q, PDF) PRINT PRINT " P1="; P PRINT " P2="; Q PRINT " Q ="; PDF GOSUB 2000 'call NPROB(X, P, Q, PDF) PRINT PRINT " P1="; P PRINT " P2="; Q PRINT " Q ="; PDF PRINT END 'of main program 500 'function alnorm(X, IUP) ' Algorithm AS66 Applied Statistics (1973) vol22 no.3 ' Evaluates the tail area of the standardised normal curve ' from x to infinity if upper is .true. or ' from minus infinity to x if upper is .false. ' Labels: 10,20,30,40 ' z,y,Temp: double ' p,q,r,a1,a2,a3,b1,b2,c1,c2,c3,c4,c5,c6: double ' d1,d2,d3,d4,d5: double ' iup:Integer (0 or 1) '*** machine dependent constants P = .398942280444#: Q = .39990348504#: r = .398942280385# A1 = 5.75885480458#: A2 = 2.62433121679#: A3 = 5.92885724438# B1 = -29.8213557807#: B2 = 48.6959930692# c1 = -3.8052E-08: c2 = .000398064794#: c3 = -.151679116635# c4 = 4.8385912808#: c5 = .742380924027#: c6 = 3.99019417011# d1 = 1.00000615302#: d2 = 1.98615381364#: d3 = 5.29330324926# d4 = -15.1508972451#: d5 = 30.789933034# z = X IF z >= zero THEN GOTO 510 IF IUP = 0 THEN IUP = 1 ELSE IUP = 0 END IF z = -z 510 IF (z <= ltone OR IUP = 1) AND z <= utzero THEN GOTO 520 Temp = zero GOTO 540 520 y = half * z * z IF z > con THEN GOTO 530 Temp = half - z * (P - Q * y / (y + A1 + B1 / (y + A2 + B2 / (y + A3)))) GOTO 540 530 TMP = (z + c1 + d1 / (z + c2 + d2 / (z + c3 + d3 / (z + c4 + d4 / (z + c5 + d5 / (z + c6)))))) Temp = r * EXP(-y) / TMP 540 IF IUP = 1 THEN alnorm = Temp ELSE alnorm = one - Temp END IF RETURN 1000 'Procedure NORMP(Z, P, Q, PDF) ' Normal distribution probabilities accurate to 1.e-15. ' Z = no. of standard deviations from the mean. ' P, Q = probabilities to the left & right of Z. P + Q = 1. ' PDF = the probability density. ' Based upon algorithm 5666 for the error function, from: ' Hart, J.F. et al, "Computer Approximations", Wiley 1968 ' Programmer: Alan Miller ' Latest revision - 30 March 1986 ' P0,P1,P2,P3,P4,P5,P6: Double ' Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7: Double ' CUTOFF,EXPNTL,ROOT2PI: Double ' ZABS:Double P0 = 220.2068679123761#: P1 = 221.2135961699311#: P2 = 112.0792914978709# P3 = 33.912866078383#: P4 = 6.37396220353165# P5 = .7003830644436881#: P6 = 3.526249659989109D-02 Q0 = 440.4137358247522#: Q1 = 793.8265125199484#: Q2 = 637.3336333788311# Q3 = 296.5642487796737#: Q4 = 86.78073220294608#: Q5 = 16.06417757920695# Q6 = 1.755667163182642#: Q7 = 8.838834764831845D-02 CUTOFF = 7.071: ROOT2PI = 2.506628274631001# z = X ZABS = ABS(z) ' |Z| > 37.0 IF ZABS > 37! THEN PDF = zero IF z > zero THEN P = one Q = zero ELSE P = zero Q = one END IF RETURN END IF ' |Z| <= 37.0 EXPNTL = EXP(-half * ZABS * ZABS) PDF = EXPNTL / ROOT2PI ' |Z| < CUTOFF = 10/sqrt(2) IF ZABS < CUTOFF THEN TMP = EXPNTL * ((((((P6 * ZABS + P5) * ZABS + P4) * ZABS + P3) * ZABS + P2) * ZABS + P1) * ZABS + P0) P = TMP / (((((((Q7 * ZABS + Q6) * ZABS + Q5) * ZABS + Q4) * ZABS + Q3) * ZABS + Q2) * ZABS + Q1) * ZABS + Q0) ' |Z| >= CUTOFF ELSE P = PDF / (ZABS + one / (ZABS + 2# / (ZABS + 3# / (ZABS + 4# / (ZABS + .65))))) END IF IF z < zero THEN Q = one - P ELSE Q = P P = one - Q END IF RETURN 2000 'Subroutine NPROB(Z, P, Q, PDF) ' P, Q = PROBABILITIES TO THE LEFT AND RIGHT OF Z ' FOR THE STANDARD NORMAL DISTRIBUTION. ' PDF = THE PROBABILITY DENSITY FUNCTION ' REFERENCE: ADAMS,A.G. AREAS UNDER THE NORMAL CURVE, ' ALGORITHM 39, COMPUTER J., VOL. 12, 197-8, 1969. ' LATEST REVISION - 23 JANUARY 1981 } ' Labels: 10,20,30 ' A0,A1,A2,A3,A4,A5,A6,A7: Double ' B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11: Double ' Y, ZABS: Double A0 = half: A1 = .398942280444#: A2 = .399903438504# A3 = 5.75885480458#: A4 = 29.8213557808# A5 = 2.62433121679#: A6 = 48.6959930692#: A7 = 5.92885724438# B0 = .398942280385#: B1 = 3.8052E-08: B2 = 1.00000615302# B3 = .000398064794#: B4 = 1.98615381364#: B5 = .151679116635# B6 = 5.29330324926#: B7 = 4.8385912808#: B8 = 15.1508972451# B9 = .742380924027#: B10 = 30.789933034#: B11 = 3.99019417011# z = X ZABS = ABS(z) IF ZABS > 12.7 THEN GOTO 2020 y = A0 * z * z PDF = EXP(-y) * B0 IF ZABS > 1.28 THEN GOTO 2010 ' Z BETWEEN -1.28 AND +1.28 Q = A0 - ZABS * (A1 - A2 * y / (y + A3 - A4 / (y + A5 + A6 / (y + A7)))) IF z < zero THEN GOTO 2030 P = one - Q RETURN ' ZABS BETWEEN 1.28 AND 12.7 2010 TMP = (ZABS - B1 + B2 / (ZABS + B3 + B4 / (ZABS - B5 + B6 / (ZABS + B7 - B8 / (ZABS + B9 + B10 / (ZABS + B11)))))) Q = PDF / TMP IF z < zero THEN GOTO 2030 P = one - Q RETURN ' Z FAR OUT IN TAIL 2020 Q = zero PDF = zero IF z < zero THEN GOTO 2030 P = one RETURN ' NEGATIVE Z, INTERCHANGE P AND Q 2030 P = Q Q = one - P RETURN 'end of file Tnormal.bas