! Last change: DOS 1 Aug 2000 8:25 pm ! *** copyright 2000 *** ! *** filename gtrou.f95 *** John F. Monahan ** ! ********************** PROGRAM PGTROU ! TEST OF ACCEPTANCE/REJECTION ALGORITHM FOR STUDENT'S T DISTRIBUTION ! let's check out 2**14 = 16384 IMPLICIT NONE INTEGER, PARAMETER :: N = 16384 INTEGER I,K,NDF REAL, DIMENSION(N) :: X REAL XS,X2S,FN,PHIXI,PHICXI,DKSP,DKSM,DKS,AD,ALPHA,DF REAL RAN,GTROU REAL(KIND=8) XI,XNMIP1,TTAIL ! 21 FORMAT(//' Sample Mean & Variance for',I12,' Student t Deviates'& & /9x,' with',I4,' degrees of freedom' & & /' 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='gtrou.out' ) ! let's check out 2**14 = 16384 FN = REAL(N) ! initialize uniform generator XS = RAN(5151917) ! do a variety of shape parameters DO K = 1,5 NDF = (1. + K*K )/2 DF = REAL(NDF) ALPHA = DF/2. ! initialize moment sums X2S = 0. XS = 0. DO I = 1,N X(I) = GTROU(DF) 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,NDF,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 XI = X(I) XNMIP1 = X(N+1-I) PHIXI = TTAIL(NDF, -XI ) PHICXI = TTAIL(NDF, XNMIP1 ) AD = AD - REAL(2*I-1)*( LOG(PHIXI) + LOG(PHICXI) )/FN - 1. 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 END DO ! LOOP ON K / ALPHA STOP END PROGRAM PGTROU REAL FUNCTION GTROU(ALPHA) ! ORIGINAL RATIO OF UNIFORMS ALGORITHM FOR GENERATING STUDENT'S T RV'S ! ! A J KINDERMAN AND J F MONAHAN (1980) "NEW METHODS FOR GENERATING ! STUDENT'S T AND GAMMA VARIABLES," COMPUTING, VOLUME 25, PP.369-377 ! IMPLICIT NONE REAL, INTENT(IN) :: ALPHA REAL U,ZZ REAL, SAVE :: ALPH = 0. REAL, SAVE :: RA,C,E,VSTAR2 REAL RAN ! SET-UP IF NECESSARY IF( ALPH .NE. ALPHA ) THEN ALPH = ALPHA ! CONSTANTS FOR QUICK ACCEPT/REJECT C = -(ALPH + 1.)/4. RA = ( (1. + 1./ALPH)**C )/4. E = 16. * RA ! CONSTANT FOR BOX VSTAR2 = 2. IF( ALPH .GT. 1. ) VSTAR2 = & & 2.*SQRT( 2.*ALPH/(ALPH-1.) ) * ( (ALPH+1.)/(ALPH-1.) )**C END IF ! ( ALPH .NE. ALPHA ) ! ! START / RESTART HERE DO ! UNRESTRICTED DO U = RAN(1) GTROU = VSTAR2*(RAN(2)-0.5)/U ZZ = GTROU*GTROU ! INNER BOUND -- QUICK ACCEPT CHECK IF( U .LE. (5.-ZZ)*RA ) RETURN ! OUTER BOUND -- QUICK REJECT CHECK ! ** BOUND NOT VALID IF ALPH < 3 ** IF( (ALPH .LT. 3 ) .OR. ( ZZ .LT. E/U-3. ) ) THEN ! FINAL RATIO OF UNIFORMS CHECK IF( U .LE. (1.+ZZ/ALPH)**C ) RETURN END IF ! (BOTH) ! REJECT -- START OVER END DO ! UNRESTRICTED DO END FUNCTION GTROU 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(KIND=8) FUNCTION TTAIL(N,T) ! COMPUTE TAIL PROBABILITIES OF STUDENT'S T DISTRIBUTION ! WITH N DEGREES OF FREEDOM PR( X > T ) ! ! N INTEGER DEGREES OF FREEDOM ** MUST BE AT LEAST 1 ** ! T REAL*8 ARGUMENT ! ! J F MONAHAN (FEB, 1990) DEPT OF STATISTICS, NCSU, RALEIGH, NC USA ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE REAL(KIND=8), INTENT(IN) :: T INTEGER, INTENT(IN) :: N REAL(KIND=8) W,B,C,SW,R,D REAL(KIND=8), PARAMETER :: PI = 3.141592653589793D0 INTEGER N2,N21,J,JM ! TTAIL = 0.5D0 IF( T .EQ. 0.D0) RETURN ! ODD OR EVEN N2 = N / 2 IF( 2*N2 .EQ. N) THEN ! DO THE EVEN CASE FIRST W = REAL(N,8)/(T*T) !KIND ! IF T FAR IN TAIL, TREAT IT CAREFULLY IF( W .GE. .625D-1 ) THEN B = 1.D0 IF( N2 .GT. 1 ) THEN C = 1.D0 N21 = N2 - 1 R = W / (1.D0 + W) DO J = 1,N21 C = C * R * REAL(2*J-1,8)/REAL(2*J,8) !KIND IF( C .LT. 1.D-15 ) EXIT B = B + C END DO ! LOOP ON J END IF ! ( N2 .EQ. 1 ) SW =SQRT( 1.D0 + W ) TTAIL = (SW-B)/(2.D0*SW) ! DONE WITH EASY EVEN CASE ELSE ! HANDLE THE CASE WHEN W IS VERY SMALL B = 0.D0 C = 1.D0 D = 1.D0 N21 = N2 - 1 R = W / (1.D0 + W) JM = MAX( 12, N21 ) ! MUST GO OUT ENOUGH TERMS TO DO SQRT(1+W) ACCURATELY DO J = 1,JM IF( J .GT. N21 ) C = 0.D0 C = C * R * REAL(2*J-1,8)/REAL(2*J,8) !KIND D = - D * W * REAL(2*J-3,8)/REAL(2*J,8) !KIND ! ADD DIFFERENCES OF LIKE POWERS B = B + (D-C) END DO ! LOOP ON J SW = SQRT( 1.D0 + W ) TTAIL = B / (2.D0*SW) ! DONE WITH HARD EVEN CASE END IF ! ( W .GT. .625D-1 ) ! DONE WITH ALL EVEN CASES ! ! TO THE CASE WHERE N IS ODD ELSE W = REAL(N,8)/(T*T) !KIND ! IF T IS FAR IN TAIL, TREAT IT CAREFULLY IF( W .GE. .625D-1 ) THEN IF( N2 .EQ. 0 ) THEN B = 0.D0 ELSE B = 1.D0 END IF ! ( N2 .EQ. 0 ) IF( N2 .GT. 1 ) THEN C = 1.D0 N21 = N2 - 1 R = W / (1.D0 + W) DO J = 1,N21 C = C * R * REAL(2*J,8)/REAL(2*J+1,8) !KIND IF( C .LT. 1.D-15 ) EXIT B = B + C END DO ! LOOP ON J END IF ! ( N2 .GT. 1 ) C = ABS(T) /SQRT(REAL(N,8)) !KIND SW = ATAN(C) TTAIL = 0.5D0 - (SW + C*B*W/(1.D0+W) )/PI ! DO SPECIAL CASE WHERE W IS SMALL ELSE B = 0.D0 IF( N .EQ. 1 ) B = 1.D0 C = 1.D0 D = -1.D0 N21 = N2 - 1 R = W / (1.D0 + W) JM = MAX( 12, N21 ) ! DO ENOUGH TERMS SO 1/2 - ATAN(X)/PI IS ACCURATE DO J = 1,JM IF( J .GT. N21 ) C = 0.D0 C = C * R * REAL(2*J,8)/REAL(2*J+1,8) !KIND D = - D * W ! TAKE DIFFERENCES OF TERMS OF LIKE POWERS B = B + (2.D0*D/REAL((2*J-1)*(2*J+1),8) - C) !KIND END DO ! LOOP ON J TTAIL = SQRT(REAL(N,8))*B/(ABS(T)*(1.D0+W)*PI) !KIND END IF ! ( W .GE. .625D-1 ) END IF ! ( N IS EVEN ) IF( T .LT. 0.D0 ) TTAIL = 1.D0 - TTAIL RETURN END FUNCTION TTAIL ! *** end of filename gtrou.f95 *********************