!********************************************************************* !* 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= 0.999999713348428 * !* P2= 2.866515718673977E-007 * !* * !* P1= 0.999999713348428 * !* P2= 2.866515718923334E-007 * !* Q= 1.486719514734297E-006 * !* * !* P1= 0.999999713348428 * !* P2= 2.866515718673977E-007 * !* Q= 1.486719514673059E-006 * !* * !* ----------------------------------------------------------------- * !* Ref.: Journal of Applied Statistics (1973) vol22 no.3. * !* * !* F90 Release By J-P Moreau, Paris. * !* (www.jpmoreau.fr) * !********************************************************************* Program TNormal real*8 P1, P2, Q, X logical UP X=5.D0 UP=.false. P1= alnorm(X,UP) UP=.true. P2= alnorm(X,UP) print *,' ' print *,' P1=', P1 print *,' P2=', P2 call NORMP(X, P1, P2, Q) print *,' ' print *,' P1=', P1 print *,' P2=', P2 print *,' Q=', Q call NPROB(X, P1, P2, Q) print *,' ' print *,' P1=', P1 print *,' P2=', P2 print *,' Q=', Q print *,' ' End ! 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. double precision function alnorm(x,upper) ! 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. double precision zero,one,half double precision con,z,y,x double precision p,q,r,a1,a2,a3,b1,b2,c1,c2,c3,c4,c5,c6 double precision d1,d2,d3,d4,d5 logical upper,up !*** machine dependent constants double precision ltone,utzero data zero/0.0d0/, one/1.0d0/, half/0.5d0/ data ltone/7.0d0/,utzero/18.66d0/ data con/1.28d0/ data p/0.398942280444d0/,q/0.39990348504d0/,r/0.398942280385d0/ data a1/5.75885480458d0/,a2/2.62433121679d0/,a3/5.92885724438d0/ data b1/-29.8213557807d0/,b2/48.6959930692d0/ data c1/-3.8052d-8/,c2/3.98064794d-4/,c3/-0.151679116635d0/ data c4/4.8385912808d0/,c5/0.742380924027d0/,c6/3.99019417011d0/ data d1/1.00000615302d0/,d2/1.98615381364d0/,d3/5.29330324926d0/ data d4/-15.1508972451d0/,d5/30.789933034d0/ up=upper z=x if(z.ge.zero)goto 10 up=.not.up z=-z 10 if(z.le.ltone.or.up.and.z.le.utzero)goto 20 alnorm=zero goto 40 20 y=half*z*z if(z.gt.con) goto 30 alnorm=half-z*(p-q*y/(y+a1+b1/(y+a2+b2/(y+a3)))) goto 40 30 alnorm=r*dexp(-y)/(z+c1+d1/(z+c2+d2/(z+c3+d3/(z+c4+d4/(z+c5+d5/(z+c6)))))) 40 if(.not.up)alnorm=one-alnorm return end SUBROUTINE 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 IMPLICIT DOUBLE PRECISION (A-H, O-Z) DATA P0, P1, P2, P3, P4, P5, P6 /220.2068679123761D0, & 221.2135961699311D0, 112.0792914978709D0, & 33.91286607838300D0, 6.373962203531650D0, & .7003830644436881D0, .3526249659989109D-01/ DATA Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7 /440.4137358247522D0, & 793.8265125199484D0, 637.3336333788311D0, & 296.5642487796737D0, 86.78073220294608D0, & 16.06417757920695D0, 1.755667163182642D0, & .8838834764831844D-1/ DATA CUTOFF /7.071D0/ DATA ROOT2PI /2.506628274631001D0/ ZABS = ABS(Z) ! |Z| > 37. IF (ZABS .GT. 37.D0) THEN PDF = 0.D0 IF (Z .GT. 0.D0) THEN P = 1.D0 Q = 0.D0 ELSE P = 0.D0 Q = 1.D0 END IF RETURN END IF ! |Z| <= 37. EXPNTL = EXP(-0.5D0*ZABS**2) PDF = EXPNTL/ROOT2PI ! |Z| < CUTOFF = 10/sqrt(2). IF (ZABS .LT. CUTOFF) THEN P = EXPNTL*((((((P6*ZABS + P5)*ZABS + P4)*ZABS + P3)*ZABS + & P2)*ZABS + P1)*ZABS + P0)/(((((((Q7*ZABS + Q6)*ZABS + & Q5)*ZABS + Q4)*ZABS + Q3)*ZABS + Q2)*ZABS + Q1)*ZABS + & Q0) ! |Z| >= CUTOFF. ELSE P = PDF/(ZABS + 1.D0/(ZABS + 2.D0/(ZABS + 3.D0/(ZABS + 4.D0/ & (ZABS + 0.65D0))))) END IF IF (Z .LT. 0.D0) THEN Q = 1.D0 - P ELSE Q = P P = 1.D0 - Q END IF RETURN END SUBROUTINE NPROB(Z,P,Q,PDF) IMPLICIT DOUBLE PRECISION (A-H,O-Z) ! 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 !******************************************************************** DATA A0,A1,A2,A3,A4,A5,A6,A7/0.5D0, 0.398942280444D0, & 0.399903438504D0, 5.75885480458D0, 29.8213557808D0, & 2.62433121679D0, 48.6959930692D0, 5.92885724438D0/, & B0,B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11/0.398942280385D0, & 3.8052D-8, 1.00000615302D0, 3.98064794D-4, 1.98615381364D0, & 0.151679116635D0, 5.29330324926D0, 4.8385912808D0, & 15.1508972451D0, 0.742380924027D0, 30.789933034D0, & 3.99019417011D0/ ZABS = ABS(Z) IF(ZABS.GT.12.7D0) GO TO 20 Y = A0*Z*Z PDF = EXP(-Y)*B0 IF(ZABS.GT.1.28D0) GO TO 10 ! Z BETWEEN -1.28 AND +1.28 Q = A0-ZABS*(A1-A2*Y/(Y+A3-A4/(Y+A5+A6/(Y+A7)))) IF(Z.LT.0.D0) GO TO 30 P = 1.D0-Q RETURN ! ZABS BETWEEN 1.28 AND 12.7 10 Q = PDF/(ZABS-B1+B2/(ZABS+B3+B4/(ZABS-B5+B6/(ZABS+B7-B8/ & (ZABS+B9+B10/(ZABS+B11)))))) IF(Z.LT.0.D0) GO TO 30 P = 1.D0-Q RETURN ! Z FAR OUT IN TAIL 20 Q = 0.D0 PDF = 0.D0 IF(Z.LT.0.D0) GO TO 30 P = 1.D0 RETURN ! NEGATIVE Z, INTERCHANGE P AND Q 30 P = Q Q = 1.D0-P RETURN END !end of file Tnormal.f90