! Last change: DOS 29 Jul 2000 12:01 pm ! *** copyright 2000 *** ! *** filename normals.f95 *** John F. Monahan ** ! ********************** PROGRAM NORMALS ! TEST PROGRAM FOR SEVERAL NORMAL ALGORITHMS ! CDFN, CDCN -- COMPUTES NORMAL CDF AND TAIL USING RATIONAL ! FUNCTION APPROXIMATIONS FROM HART BOOK ! ALNPHI -- COMPUTES NATURAL LOG OF NORMAL CDF USING ! RATIONAL FUNCTION APPROX (MONAHAN, 1981) ! CDFVED,CDCVED -- COMPUTES NORMAL CDF AND TAIL USING A SIMPLE ! APPROX BY VEDDER -- QUICK AND DIRTY ! CDFN, CDCN ARE MOST ACCURATE ! ALNPHI IS NEXT, VERY CLOSE, BUT DESIGNED FOR DIFFERENT JOB ! CDFVED,CDCVED ARE SIMPLE, GOOD FOR A CALCULATOR, BUT NOT LIBRARY ! IMPLICIT NONE REAL X,EF,EC,EFM,ECM,EFV,ECV REAL AL10,CDFVED,CDCVED,ALNPHI,CDFN,CDCN,RAN INTEGER I,J ! 21 FORMAT(F8.3,3F11.7,5X,3F11.7) 22 FORMAT(I4,2F11.5) 25 FORMAT('FIRST TEST -- COMPARE TO TABLE VALUES') 26 FORMAT(' X E(ALNPHI(X)) CDFN(X) CDFVED(X)',5X, & & 'E(ALNPHI(-X)) CDCN(X) CDCVED(X)') 27 FORMAT(//'JUST CHECK THE TAIL'/' X ALNPHI(-X)') 28 FORMAT(//'LAST TEST -- RANDOM POINTS -- COMPARE TO EACH OTHER') ! OPEN ( UNIT=6, FILE='normals.out' ) ! WRITE(6,25) WRITE(6,26) ! FIRST REPRODUCE A TABLE DO I = 1,28 X = REAL(I-1)/5. EFV = CDFVED(X) ECV = CDCVED(X) EF = EXP( ALNPHI( X ) ) EC = EXP( ALNPHI( -X ) ) EFM = CDFN(X) ECM = CDCN(X) WRITE(6,21) X,EF,EFM,EFV,EC,ECM,ECV END DO ! LOOP ON I ! ! NEXT A TAIL TABLE, SEE A&S TABLE 26.2 WRITE(6,27) AL10 = LOG(10.) DO I = 5,100,5 X = REAL(I) EF = - ALNPHI(-X)/AL10 WRITE(6,22) I,EF END DO ! LOOP ON I ! ! LAST JUST SOME RANDOM VALUES -- TEST BY COMPARING ! ! INITIALIZE SEED FOR UNIFORM GENERATOR X = RAN(5151917) ! WRITE(6,28) WRITE(6,26) ! LAST JUST SOME RANDOM VALUES -- TEST BY COMPARING DO J = 1,5 DO I = 1,5 X = REAL(J)*4.*(RAN(I)-.5) EFV = CDFVED(X) ECV = CDCVED(X) EF = EXP( ALNPHI( X ) ) EC = EXP( ALNPHI( -X ) ) EFM = CDFN(X) ECM = CDCN(X) WRITE(6,21) X,EF,EFM,EFV,EC,ECM,ECV END DO ! LOOP ON I END DO ! LOOP ON J STOP END PROGRAM NORMALS REAL FUNCTION CDFVED(X) ! NORMAL DISTRIBUTION FUNCTION USING A SIMPLE APPROXIMATION ! BY J.D. VEDDER (1993) COMP STAT & DATA ANAL PP.119-123 ! THIS VERSION IS THE CDF IMPLICIT NONE REAL, INTENT(IN) :: X REAL, PARAMETER :: A = 1.595769 ! SQRT(8/PI) REAL, PARAMETER :: B = .0726712 ! SQRT(2/PI)*(4-PI)/(3*PI) REAL XX,U XX = X*X U = X * (A + B*XX) ! APPROXIMATION IS MODIFICATION OF LOGISTIC CDFVED = 1./(1. + EXP(-U) ) RETURN END FUNCTION CDFVED REAL FUNCTION CDCVED(X) ! NORMAL DISTRIBUTION FUNCTION USING A SIMPLE APPROXIMATION ! BY J.D. VEDDER (1993) COMP STAT & DATA ANAL PP.119-123 ! THIS VERSION IS THE COMPLEMENT OR TAIL PROB IMPLICIT NONE REAL, INTENT(IN) :: X REAL, PARAMETER :: A = 1.595769 ! SQRT(8/PI) REAL, PARAMETER :: B = .0726712 ! SQRT(2/PI)*(4-PI)/(3*PI) REAL XX,U XX = X*X U = X * (A + B*XX) ! APPROXIMATION IS MODIFICATION OF LOGISTIC CDCVED = 1./(1. + EXP(U) ) RETURN END FUNCTION CDCVED REAL FUNCTION CDFN(X) ! CDF FOR NORMAL DIST 22 MAY 74 JFM ! ** RECODED AS SEPARATE SUBPROGRAMS ** 14 JULY 1988 ** ! ** RECODED FOR READABILITY ** 25 JUNE 1996 ** ! ** RECODED FOR FORTRAN 95 ** 30 JAN 2000 ** ! ! CDFN=CDF FOR NORMAL DIST CDCN=1-CDFN ! CDFN=0.5*(1.0+ERF()/SQRT(2))) ! ! RANGE FOR X METHOD SIZE NO. ACCURACY ! (0.0,0.14) CONTINUED FRACTION 3 6.7.9 7 ! (0.14,3.54) RATIONAL FUNCTION 3,4 5663 6,MOSTLY 7 ! (3.54,11.3) RATIONAL FUNCTION 4,5 5704 5,MOSTLY 6 ! (11.3,13.15) RATIONAL FUNCTION 1,2 5721 5 ! NO. REFERS TO EQN NO. OR METHOD FROM 'COMPUTER APPROXIMATIONS' ! ED BY J F HART WILEY 1968 ! REFERENCE FOR ACCURACY 'TABLES OF ERROR FN AND ITS DERIVATIVE' ! NATL BUREAU OF STANDARDS AMS 41 1958 ! ACCURACY IS IN SIG FIGS FOR UNIVAC 1108 ! IMPLICIT NONE REAL, INTENT(IN) :: X REAL, DIMENSION(4) :: P2 REAL, DIMENSION(5) :: Q2,P3 REAL, DIMENSION(6) :: Q3 REAL, DIMENSION(2) :: P4 REAL, DIMENSION(3) :: Q4 REAL XA,XS,V DATA P2/ .100046412E+2, .842655286E+1, .346025933E+1, & & .562353612/ DATA Q2/ .100046412E+2, .197155807E+2, .157022881E+2, & & .609074879E+1, 1.0 / DATA P3/ .183983860E+2, .224003952E+2, .130613857E+2, & & .402838962E+1, .564201377/ DATA Q3/ .183983860E+2, .431607573E+2, .433645871E+2, & & .236403038E+2, .714083720E+1, 1.0/ DATA P4/ .148459208, .564187743/ DATA Q4/ .511437285, .262770594, 1.0/ ! DEFINE SOME STATEMENT FUNCTIONS ! RESCALE SINCE ORIGINAL CODING WAS FOR ERROR FUNCTION XA = ABS(X/1.414214) ! ABS( X / SQRT(2) ) XS = X*X/2. IF( XA .LE. 0.1 ) THEN ! FIRST PIECE FOR SMALLEST X -- USE CONTINUED FRACTION V = 1.12837917*XA / ( EXP(XS) * (1.0-2.0*XS/(3.0+0.8*XS)) ) CDFN = ( 1. + SIGN( V, X) )/2. RETURN END IF ! IF( XA .LE. 0.1 ) IF( XA .LE. 2.5 ) THEN ! SECOND PIECE -- USE STRAIGHT RATIONAL FUNCTION (3,4) CDFN = (P2(1)+XA*(P2(2)+XA*(P2(3)+XA*P2(4)))) / & & (2.*EXP(XS)*(Q2(1)+XA*(Q2(2)+XA*(Q2(3)+XA*(Q2(4)+XA*Q2(5)))))) IF( X .GT. 0. ) CDFN = 1. - CDFN RETURN END IF ! IF( XA .LE. 2.5 ) IF( X .GT. 5. ) THEN ! TOO FAR IN RIGHT TAIL -- RETURN ONE CDFN = 1. RETURN END IF ! IF( X .GT. 5. ) IF( XA .LE. 8. ) THEN ! THIRD PIECE -- RATIONAL FUNCTION (4,5) FOR NEAR TAIL CDFN = (P3(1)+XA*(P3(2)+XA*(P3(3)+XA*(P3(4)+XA*P3(5))))) / (2.* & & (Q3(1)+XA*(Q3(2)+XA*(Q3(3)+XA*(Q3(4)+XA*(Q3(5)+XA*Q3(6))))))) CDFN = CDFN / EXP(XS) IF( X .GT. 0. ) CDFN = 1. - CDFN RETURN END IF ! IF( XA .LE. 8. ) IF( X .GE. -9.3 ) THEN ! FOURTH (TAIL) PIECE -- ANOTHER RATIONAL FUNCTION (1,2) FOR TAIL CDFN = (P4(1)+XA*P4(2))*EXP(-XS)/(2.*(Q4(1)+XA*(Q4(2)+XA*Q4(3)))) RETURN END IF ! IF( XA .LE. 9.3 ) ! TOO FAR IN LEFT TAIL -- RETURN ZERO CDFN = 0. RETURN END FUNCTION CDFN REAL FUNCTION CDCN(X) ! 1 - CDF FOR NORMAL DIST 22 MAY 74 JFM ! ** RECODED AS SEPARATE SUBPROGRAMS ** 14 JULY 1988 ** ! ** RECODED FOR READABILITY ** 25 JUNE 1996 ** ! ** RECODED FOR FORTRAN 95 ** 30 JAN 2000 ** ! ! CDFN=CDF FOR NORMAL DIST CDCN=1-CDFN ! CDFN=0.5*(1.0+ERF()/SQRT(2))) ! ! RANGE FOR X METHOD SIZE NO. ACCURACY ! (0.0,0.14) CONTINUED FRACTION 3 6.7.9 7 ! (0.14,3.54) RATIONAL FUNCTION 3,4 5663 6,MOSTLY 7 ! (3.54,11.3) RATIONAL FUNCTION 4,5 5704 5,MOSTLY 6 ! (11.3,13.15) RATIONAL FUNCTION 1,2 5721 5 ! NO. REFERS TO EQN NO. OR METHOD FROM 'COMPUTER APPROXIMATIONS' ! ED BY J F HART WILEY 1968 ! REFERENCE FOR ACCURACY 'TABLES OF ERROR FN AND ITS DERIVATIVE' ! NATL BUREAU OF STANDARDS AMS 41 1958 ! ACCURACY IS IN SIG FIGS FOR UNIVAC 1108 ! IMPLICIT NONE REAL, INTENT(IN) :: X REAL, DIMENSION(4) :: P2 REAL, DIMENSION(5) :: Q2,P3 REAL, DIMENSION(6) :: Q3 REAL, DIMENSION(2) :: P4 REAL, DIMENSION(3) :: Q4 REAL XA,XS,V DATA P2/ .100046412E+2, .842655286E+1, .346025933E+1, & & .562353612/ DATA Q2/ .100046412E+2, .197155807E+2, .157022881E+2, & & .609074879E+1, 1.0 / DATA P3/ .183983860E+2, .224003952E+2, .130613857E+2, & & .402838962E+1, .564201377/ DATA Q3/ .183983860E+2, .431607573E+2, .433645871E+2, & & .236403038E+2, .714083720E+1, 1.0/ DATA P4/ .148459208, .564187743/ DATA Q4/ .511437285, .262770594, 1.0/ ! DEFINE SOME STATEMENT FUNCTIONS ! RESCALE SINCE ORIGINAL CODING WAS FOR ERROR FUNCTION XA = ABS(X/1.414214) ! ABS( X / SQRT(2) ) XS = X*X/2. IF( XA .LE. 0.1 ) THEN ! FIRST PIECE FOR SMALLEST X -- USE CONTINUED FRACTION V = 1.12837917*XA / ( EXP(XS) * (1.0-2.0*XS/(3.0+0.8*XS)) ) CDCN = ( 1. - SIGN( V, X) )/2. RETURN END IF ! IF( XA .LE. 0.1 ) IF( XA .LE. 2.5 ) THEN ! SECOND PIECE -- USE STRAIGHT RATIONAL FUNCTION (3,4) CDCN = (P2(1)+XA*(P2(2)+XA*(P2(3)+XA*P2(4)))) / & & (2.*EXP(XS)*(Q2(1)+XA*(Q2(2)+XA*(Q2(3)+XA*(Q2(4)+XA*Q2(5)))))) IF( X .LT. 0. ) CDCN = 1. - CDCN RETURN END IF ! IF( XA .LE. 2.5 ) IF( X .LT. -5. ) THEN ! TOO FAR IN LEFT TAIL -- RETURN ONE CDCN = 1. RETURN END IF ! IF( X .GT. 5. ) IF( XA .LE. 8. ) THEN ! THIRD PIECE -- RATIONAL FUNCTION (4,5) FOR NEAR TAIL CDCN = (P3(1)+XA*(P3(2)+XA*(P3(3)+XA*(P3(4)+XA*P3(5))))) / (2.* & & (Q3(1)+XA*(Q3(2)+XA*(Q3(3)+XA*(Q3(4)+XA*(Q3(5)+XA*Q3(6))))))) CDCN = CDCN / EXP(XS) IF( X .LT. 0. ) CDCN = 1. - CDCN RETURN END IF ! IF( XA .LE. 8. ) IF( X .LE. 9.3 ) THEN ! FOURTH (TAIL) PIECE -- ANOTHER RATIONAL FUNCTION (1,2) FOR TAIL CDCN = (P4(1)+XA*P4(2))*EXP(-XS)/(2.*(Q4(1)+XA*(Q4(2)+XA*Q4(3)))) RETURN END IF ! IF( XA .LE. 9.3 ) ! TOO FAR IN RIGHT TAIL -- RETURN ZERO CDCN = 0. RETURN END FUNCTION CDCN REAL FUNCTION ALNPHI(X) ! COMPUTES THE NATURAL LOGARITHM OF THE NORMAL DF ! J F MONAHAN JAN 1981 DEPT. OF STAT., N C S U, RALEIGH, N C 27650 ! RELATIVE ACCURACY 4E-6 FOR -4 14 IF( X .GT. 14.) RETURN AX = ABS(X) IF( AX .GT. 2. ) THEN ! DO THE TAILS FIRST XX = X*X XR32 = (((PR3/XX+PR2)/XX+PR1)/AX+AX)/((QR2/XX+QR1)/XX+1.) IF( X .GT. 0. ) THEN OMPHI = EXP(-X*X/2.)*RSR2PI/XR32 ! OMPHI=1-PHI NEXT USE PADE APPROX TO GET LN(1-OMPHI) ALNPHI = -OMPHI*(6-OMPHI)/(6-4*OMPHI) RETURN END IF ! IF( X .GT. 0. ) ALNPHI = LOG(RSR2PI/XR32) - X*(0.5*X) RETURN END IF ! IF( AX .GT. 2. ) IF( X .LE. 0.60 ) THEN ! NOW DO THE RANGE -2 < X < 0.6 ALNPHI = (PS0+X*(PS1+X*(PS2+X*(PS3+X*PS4))))/(1.0+X*(QS1+X*QS2)) RETURN END IF ! IF( X .LT. 0.60 ) ! LAST TAKE CARE OF 0.6 < X < 2 ALNPHI = (PT0+X*(PT1+X*(PT2+X*(PT3+X*PT4)))) & & *EXP(-0.5*X*X)/(1.0+X*(QT1+X*QT2)) RETURN END FUNCTION ALNPHI REAL FUNCTION RAN(IDUM) ! PORTABLE IMPLEMENTATION OF UNIFORM PSEUDORANDOM NUMBER GENERATOR ! LEWIS, GOODMAN, & MILLER MULTIPLICATIVE GENERATOR ! X(N+1) = MOD( 16807*X(N), 2**31-1 ) ! ! P. BRANTLEY, B.L. FOX, L. SCHRAGE (1983) A GUIDE TO SIMULATION ! SPRINGER-VERLAG, NEW YORK. PP. 200-202. ! UPDATED VERSION OF ! LINUS SCHRAGE (1979)'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR' ! ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 5, PP. 132-138. ! ! ARGUMENT ! IDUM INTEGER FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS ! IMPLICIT NONE INTEGER, INTENT(IN) :: IDUM REAL, PARAMETER :: RP = 4.656612875E-10 ! 1/P INTEGER, PARAMETER :: A = 16807 ! MULTIPLIER INTEGER, PARAMETER :: B = 127773 ! B = P / A INTEGER, PARAMETER :: C = 2836 ! C = P MOD A INTEGER, PARAMETER :: P = 2147483647 ! MODULUS 2**31 - 1 INTEGER, SAVE :: IX = 0 INTEGER K1 ! ! IF NOT FIRST CALL, THEN SKIP SETTING SEED IF( IX .EQ. 0) IX = IDUM ! WRITE NUMBER AS ALPHA*2**16 + BETA K1 = IX / B IX = A*( IX - K1*B) - K1*C ! ABOVE DOES A*IX MOD B -K1*C IF( IX .LT. 0 ) IX = IX + P ! RP BELOW IS RECIPROCAL OF P RAN = REAL(IX)*RP RETURN END FUNCTION RAN ! *** end of filename normals.f95 *********************