! Last change: DOS 29 Jul 2000 12:28 pm ! *** copyright 2000 *** ! *** filename ttail.f95 *** John F. Monahan ** ! ********************** ! ! SOME COMPILERS REQUIRE THAT MODULES ARE AVAILABLE BEFORE ! REQUESTED, OR AT THE BEGINNING ! MODULE PASSINFO ! N AND P GIVEN VALUES IN PROGRAM REAL(KIND=8) :: P ! AND ACCESSED BY EXTERNAL FUNCTION G INTEGER :: N ! WHICH IS CALLED BY SUBPROGRAM SECANT END MODULE PASSINFO ! WHICH IS UNAWARE OF MODULE ! PROGRAM PTTAIL ! TEST OF CODE TO COMPUTE TAIL PROBABILITIES OF STUDENT'S T ! COMPUTE PERCENTILE POINTS ACCORDING TO TABLE 26.10 IN A & S USE PASSINFO IMPLICIT NONE REAL(KIND=8), EXTERNAL :: G REAL(KIND=8) TTAIL REAL(KIND=8), DIMENSION(7) :: PTABLE,XTABLE REAL(KIND=8) X1,X2,PC INTEGER J ! COMPUTE CRITICAL VALUES FOR THESE P'S DATA PTABLE/ .4D0, .25D0, .1D0, .025D0, .01D0, .001D0, .5D-5 / ! 20 FORMAT(' TABLE OF T CRITICAL VALUES -- ONE-SIDED PROBS ARE'& & /' N',4X,F4.2,4F11.3,F11.4,F11.6/) 21 FORMAT(I4,4F15.10) 22 FORMAT(I3,F10.6,6F11.6) ! OPEN OUTPUT FILE OPEN( UNIT=6, FILE='ttail.out' ) WRITE(6,20) PTABLE ! FOR DIFFERENT VALUES OF DF N DO N = 1,60 IF( ( N .LT. 10 ) .OR. ( MOD(N,5) .EQ. 0 ) ) THEN ! FOR DIFFERENT VALUES OF P DO J = 1,7 P = PTABLE(J) ! COMMON STARTING VALUES -- HOPE THEY'RE GOOD ENOUGH X1 = .5D0 X2 = 3.5D0 ! FIND CRITICAL VALUE CALL SECANT(G,X1,X2,30,1.D-9,1.D-9) ! CHECK IT ANYWAY PC = TTAIL(N,X1) IF( ABS( (PC-P)/P ) .GT. 1.D-8 ) WRITE(6,21) N,X1,X2,P,PC XTABLE(J) = X1 ! SAVE CRITICAL VALUE IN TABLE END DO ! LOOP ON J WRITE(6,22) N,XTABLE END IF ! END DO ! LOOP ON N STOP END PROGRAM PTTAIL REAL(KIND=8) FUNCTION G(X) USE PASSINFO IMPLICIT NONE REAL(KIND=8), INTENT(IN) :: X REAL(KIND=8) TTAIL ! RELATIVE ERROR TO GET THE FAR TAIL G = ( TTAIL(N,X) - P )/P RETURN END FUNCTION G 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 SUBROUTINE SECANT(F,X1,X2,MIT,EPSX,EPSF) ! THE SECANT METHOD FOR FINDING ROOTS TO F ! RETURNS ROOT IN X1, VALUE OF F THERE IN X2 ! EPSF: STOP IF ABS(F(X)) .LE. EPSF ! EPSX: STOP IF CHANGE IN X IS .LE. EPSX ! ! *** DOUBLE PRECISION VERSION *** IMPLICIT NONE REAL(KIND=8) F REAL(KIND=8), INTENT(IN OUT) :: X1,X2 REAL(KIND=8), INTENT(IN) :: EPSX,EPSF INTEGER, INTENT(IN) :: MIT REAL(KIND=8) ESPX,ESPF,FX1,FX2,XN,FXN INTEGER IT ! DON'T LET EPSILONS BE TOO SMALL ESPX=MAX(1.D-14,EPSX) ESPF=MAX(1.D-14,EPSF) ! TWO STARTING POINTS FX1=F(X1) FX2=F(X2) DO IT = 1,MIT ! FIND NEW POINT XN=X2-FX2*(X2-X1)/(FX2-FX1) FXN=F(XN) ! TEST FOR CONVERGENCE IF( ( ABS(FXN) .LT. ESPF) .OR. (ABS(XN-X2) .LT. ESPX) ) EXIT ! UPDATE X1=X2 FX1=FX2 X2=XN FX2=FXN END DO ! LOOP ON IT X1=XN X2=FXN RETURN END SUBROUTINE SECANT ! *** end of filename ttail.f95 *********************