! Last change: DOS 1 Aug 2000 8:15 pm ! *** copyright 2000 *** ! *** filename gnpolr.f95 *** John F. Monahan ** ! ********************** PROGRAM PGNPOLR ! TEST OF POLAR METHOD FOR GENERATING NORMAL RANDOM VARIABLES ! let's check out 2**14 = 16384 IMPLICIT NONE INTEGER, PARAMETER :: N = 16384 INTEGER I REAL, DIMENSION(N) :: X REAL XS,X2S,FN,PHIXI,DKSP,DKSM,DKS,AD REAL RAN,GNPOLR,CDFN,ALNPHI ! 21 FORMAT(' Sample Mean and Variance for ',I8,' Normal Deviates' & & /' Mean is',F16.8,4x,'Variance is',F16.8) 22 FORMAT(' Goodness of Fit Statistics: '/4x,'Anderson-Darling is',& &F12.6/4x,'Compare to .05, .01 critical values of 2.492, 3.857'//& & 4x,'Kolmogorov-Smirnov D+',F9.6,4x,' D-',F9.6/ & & 4x,'SQRT(N)*MAX(D+,D-) is',F9.4/4x,'Compare to .05, .01', & & ' critical values of 1.36, 1.63') ! ! open output file OPEN( UNIT=6, FILE='gnpolr.out' ) ! let's check out 2**14 = 16384 FN = REAL(N) ! initialize uniform generator XS = RAN(5151917) ! initialize moment sums X2S = 0. XS = 0. DO I = 1,N X(I) = GNPOLR(I) XS = XS + X(I) IF( I .GT. 1 ) X2S = X2S + (XS - I*X(I))*(XS - I*X(I))/(I*(I-1)) END DO ! LOOP ON I ! get sample mean and variance XS = XS/FN X2S = X2S/(FN-1.) WRITE(6,21) N,XS,X2S ! sort the data -- get order statistics CALL HSORT(X,N) ! now compute K-S and A-D statistics DKSP = 0. DKSM = 0. AD = 0. DO I = 1,N AD = AD - (2*I-1)*( ALNPHI(X(I)) + ALNPHI(-X(N+1-I)) )/FN - 1. PHIXI = CDFN( X(I) ) DKSP = MAX( DKSP, REAL(I)/FN - PHIXI ) DKSM = MAX( DKSM, PHIXI - REAL(I-1)/FN ) END DO ! LOOP ON I ! usual standardization of K-S DKS = SQRT(FN)*MAX( DKSP, DKSM ) WRITE(6,22) AD,DKSP,DKSM,DKS STOP END PROGRAM PGNPOLR REAL FUNCTION GNPOLR(IXX) ! MARSAGLIA AND BRAY POLAR METHOD FOR GENERATING NORMAL(0,1) RV'S ! ! G. MARSAGLIA & T. A. BRAY (1964) "A CONVENIENT METHOD FOR GENERATING ! NORMAL VARIABLES," SIAM REVIEW, VOLUME 6, PP. 260-264. ! ! ARGUMENT IS DUMMY -- NEVER USED ! PAIRS OF NORMALS ARE PRODUCED ON ODD-NUMBERED CALLS, EVEN-NUMBERED ! CALLS JUST RETURN THE SECOND ONE ! IMPLICIT NONE INTEGER, INTENT(IN) :: IXX REAL, SAVE :: U2 REAL, PARAMETER :: TWOPI = 6.2831853 REAL U1,W REAL RAN INTEGER, SAVE :: IND = 1 ! BEFORE FIRST CALL, INDICATOR SAYS LAST (0) WAS EVEN ! CHANGE ODD/EVEN IND = -IND ! IF( IND .GT. 0 ) THEN ! EVEN-NUMBERED CALL -- TAKE ONE THAT'S WAITING GNPOLR = U2 ELSE ! ODD CALLS START/RESTART HERE DO ! NOTE UNRESTRICTED DO U1 = RAN(1) U2 = RAN(2) U1 = U1 + U1 - 1. U2 = U2 + U2 - 1. ! (U1, U2) POINT UNIFORM IN BOX W = U1*U1 + U2*U2 ! TEST TO SEE IF POINT IS IN THE CIRCLE IF( W .LE. 1. ) EXIT END DO ! EXIT (ACCEPT) LEAVES UNRESTRICTED DO ! W IS RADIUS W = SQRT( -2.*LOG(W)/W ) ! SAVE U2 FOR NEXT CALL AND RETURN U1*W U2 = U2*W GNPOLR = U1*W END IF ! ( IND .GT. 0 ) RETURN END FUNCTION GNPOLR SUBROUTINE HSORT(K,N) ! HEAPSORT ALGORITHM FOR SORTING ON VECTOR OF KEYS K OF LENGTH N ! ! ARGUMENTS ! K REAL VECTOR OF KEYS TO BE SORTED ! N NUMBER OF ITEMS TO BE SORTED ! ! TO SORT A PARALLEL VECTOR OF RECORDS, USE HKSORT ! ! J F MONAHAN (DEC, 1999) FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION(N), INTENT(IN OUT) :: K REAL KK INTEGER I,L,NCUR ! ! DO NOTHING IF THERE'S NOTHING TO DO IF( N .LE. 1 ) RETURN ! INITIALIZE TO BUILDHEAP PART (LOOP ON L) L = N/2 + 1 NCUR = N DO I = L,1,-1 CALL HEAPIFY(I) END DO ! LOOP ON I DO I = 2,N ! SWITCH CURRENT LARGEST WITH BOTTOM KK = K(1) K(1) = K(NCUR) K(NCUR) = KK ! REHEAP WITH ONE SHORTER NCUR = NCUR - 1 CALL HEAPIFY(1) END DO ! LOOP ON NCUR RETURN CONTAINS SUBROUTINE HEAPIFY(II) INTEGER, INTENT(IN) :: II INTEGER I,J I = II DO J = 2*I ! IS IT A LEAF OR ARE THERE SONS? IF( J > NCUR ) EXIT ! A LEAF IF( J < NCUR ) THEN ! ANOTHER SON OF I IF( K(J+1) > K(J) ) J = J+1 ! LARGER SON IS K END IF ! A LEAF IF( K(J) > K(I) ) THEN ! EXCHANGE KK = K(J) K(J) = K(I) K(I) = KK I = J ELSE ! EXIT -- HEAP PROPERTY EXIT END IF END DO ! WHILE NOT A LEAF END SUBROUTINE HEAPIFY END SUBROUTINE HSORT 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 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 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 ! *** end of filename gnpolr.f95 *********************