! Last change: DOS 2 Aug 2000 9:24 pm ! *** copyright 2000 *** ! *** filename perfect.f95 *** John F. Monahan ** ! ********************** PROGRAM PERFECT ! TEST OF BOX-MULLER METHOD FOR GENERATING NORMAL RANDOM VARIABLES ! *** BUT USE HALTON SCHEME TO GENERATE "TOO PERFECT" SAMPLE *** ! ! 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,GNBXML,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='perfect.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) = GNBXML(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 PERFECT REAL FUNCTION RAN(I) ! *** NOT THE USUAL UNIFORM RANDOM NUMBER GENERATOR *** ! GENERATE TWO DIMENSIONAL QUASIRANDOM POINTS UNIFORMLY ! DISTRIBUTED IN (0,1)x(0,1) USING HALTON SCHEME ! IMPLICIT NONE INTEGER, INTENT(IN) :: I INTEGER, PARAMETER :: BASE2 = 2 ! BASE FOR ONE SEQUENCE INTEGER, PARAMETER :: BASE3 = 3 ! BASE FOR OTHER SEQUENCE INTEGER, SAVE :: NNZ2 = 1 INTEGER, SAVE :: NNZ3 = 1 INTEGER, SAVE, DIMENSION(100) :: ILIST2 = 0 INTEGER, SAVE, DIMENSION(100) :: ILIST3 = 0 ! ! INTERFACE BLOCK INTERFACE subroutine ncrmnt(list,base,nnzero) integer, dimension(:) :: list integer, INTENT(IN) :: base integer, INTENT(IN OUT) :: nnzero end subroutine ncrmnt real function expand(list,base,nnzero) integer, dimension(:) :: list integer, INTENT(IN) :: base integer, INTENT(IN OUT) :: nnzero end function expand END INTERFACE ! ! HALTON SCHEME VALUES, USE BASES 2 AND 3 ! IF( I .EQ. 1 ) THEN ! IF I = 1 THEN GET NEXT POINT CALL NCRMNT( ILIST2, BASE2, NNZ2) RAN = EXPAND( ILIST2, BASE2, NNZ2) ELSE CALL NCRMNT( ILIST3, BASE3, NNZ3) RAN = EXPAND( ILIST3, BASE3, NNZ3) END IF ! ( I .EQ. 1 ) RETURN END FUNCTION RAN REAL FUNCTION GNBXML(IXX) ! BOX-MULLER METHOD FOR GENERATING NORMAL(0,1) RV'S ! ! GEORGE E. P. BOX & MERVIN MULLER (1958) "A NOTE ON THE GENERATION ! RANDOM NORMAL DEVIATES," ANNALS OF MATHEMATICAL STATISTICS, ! VOLUME 29, PP. 610-611. ! ! 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 :: U REAL, PARAMETER :: TWOPI = 6.2831853 REAL V,W REAL RAN INTEGER, SAVE :: IND = 1 ! BEFORE FIRST CALL, INDICATOR SAYS LAST (0) WAS EVEN ! CHANGE ODD/EVEN IND = -IND ! IF EVEN (IND IS POSITIVE) THEN GRAB WAITING ONE IF( IND .GT. 0 ) THEN ! EVEN-NUMBERED CALL -- TAKE ONE THAT'S WAITING GNBXML = U ELSE ! ODD CALLS START/RESTART HERE U = RAN(1) V = RAN(2) ! W IS RADIUS -- CHI WITH 2 DF W = SQRT( -2.*LOG(V) ) GNBXML = SIN( TWOPI*U ) * W U = COS( TWOPI*U ) * W ! SAVE ONE FOR NEXT CALL AND RETURN FIRST ONE END IF ! ( IND .GT. 0 ) END IF RETURN END FUNCTION GNBXML subroutine ncrmnt(list,base,nnzero) ! increments by one a number expressed as an integer in base "base" ! with each digit residing in "list" ! arguments ! list integer vector of length at least nnzero holding digit ! of base "base" expansion -- dummy dimensioned as length 1 ! base integer base ! nnzero integer giving the number of nonzero digits ! ! if list initialized to all zeros, then after n calls to nrcrmnt ! nnzero ! n = SUM list(j)*base**(j) ! j=1 ! implicit none integer, dimension(:) :: list integer, INTENT(IN) :: base integer, INTENT(IN OUT) :: nnzero integer j ! if nnzero is goofy, then do nothing if( nnzero .lt. 0 ) return ! increment first digit j = 1 do ! repeat loop without condition if( j .gt. nnzero ) nnzero = j list(j) = list(j) + 1 ! if no carry, then done if( list(j) .lt. base ) exit ! carrying list(j) = 0 j = j + 1 end do ! end loop return end subroutine ncrmnt real function expand(list,base,nnzero) ! computes the real number corresponding to its base "base" ! expansion to nnzero digits listed in list ! ! nnzero ! x = expand(list,base,nnzero) = SUM list(j)*base**(-j) ! j=1 ! arguments ! list integer vector of length at least nnzero holding digit ! of base "base" expansion -- dummy dimensioned as length 1 ! base integer base ! nnzero integer giving the number of nonzero digits ! ! ** note that all digits are to the right of the radix point ** ! ** so that expand(list,base,nnzero) is between 0 and 1 ** ! implicit none integer, dimension(:) :: list integer, INTENT(IN) :: base integer, INTENT(IN OUT) :: nnzero integer j,jj real flbase ! get floating point value of base flbase = real(base) ! initialize value to zero expand = 0. ! increment j through number of digits given do j = 1,nnzero ! do the least significant first jj = nnzero + 1 - j expand = list(jj) + expand/flbase end do ! loop on j ! last divide expand = expand / base return end function expand 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 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 perfect.f95 *********************