SUBROUTINE TPPF(P,NU,PPF) C C PURPOSE--THIS SUBROUTINE COMPUTES THE PERCENT POINT C FUNCTION VALUE FOR THE STUDENT'S T DISTRIBUTION C WITH INTEGER DEGREES OF FREEDOM PARAMETER = NU. C THE STUDENT'S T DISTRIBUTION USED C HEREIN IS DEFINED FOR ALL X, C AND ITS PROBABILITY DENSITY FUNCTION IS GIVEN C IN THE REFERENCES BELOW. C NOTE THAT THE PERCENT POINT FUNCTION OF A DISTRIBUTION C IS IDENTICALLY THE SAME AS THE INVERSE CUMULATIVE C DISTRIBUTION FUNCTION OF THE DISTRIBUTION. C INPUT ARGUMENTS--P = THE SINGLE PRECISION VALUE C (BETWEEN 0.0 (EXCLUSIVELY) C AND 1.0 (EXCLUSIVELY)) C AT WHICH THE PERCENT POINT C FUNCTION IS TO BE EVALUATED. C --NU = THE INTEGER NUMBER OF DEGREES C OF FREEDOM. C NU SHOULD BE POSITIVE. C OUTPUT ARGUMENTS--PPF = THE SINGLE PRECISION PERCENT C POINT FUNCTION VALUE. C OUTPUT--THE SINGLE PRECISION PERCENT POINT FUNCTION . C VALUE PPF FOR THE STUDENT'S T DISTRIBUTION C WITH DEGREES OF FREEDOM PARAMETER = NU. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--NU SHOULD BE A POSITIVE INTEGER VARIABLE. C --P SHOULD BE BETWEEN 0.0 (EXCLUSIVELY) C AND 1.0 (EXCLUSIVELY). C OTHER DATAPAC SUBROUTINES NEEDED--NORPPF. C FORTRAN LIBRARY SUBROUTINES NEEDED--DSIN, DCOS, DSQRT, DATAN. C MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION. C LANGUAGE--ANSI FORTRAN. C COMMENT--FOR NU = 1 AND NU = 2, THE PERCENT POINT FUNCTION C FOR THE T DISTRIBUTION EXISTS IN SIMPLE CLOSED FORM C AND SO THE COMPUTED PERCENT POINTS ARE EXACT. C --FOR OTHER SMALL VALUES OF NU (NU BETWEEN 3 AND 6, C INCLUSIVELY), THE APPROXIMATION C OF THE T PERCENT POINT BY THE FORMULA C GIVEN IN THE REFERENCE BELOW IS AUGMENTED C BY 3 ITERATIONS OF NEWTON'S METHOD FOR C ROOT DETERMINATION. C THIS IMPROVES THE ACCURACY--ESPECIALLY FOR C VALUES OF P NEAR 0 OR 1. C REFERENCES--NATIONAL BUREAU OF STANDARDS APPLIED MATHMATICS C SERIES 55, 1964, PAGE 949, FORMULA 26.7.5. C --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE C DISTRIBUTIONS--2, 1970, PAGE 102, C FORMULA 11. C --FEDERIGHI, 'EXTENDED TABLES OF THE C PERCENTAGE POINTS OF STUDENT'S T C DISTRIBUTION, JOURNAL OF THE C AMERICAN STATISTICAL ASSOCIATION, C 1969, PAGES 683-688. C --HASTINGS AND PEACOCK, STATISTICAL C DISTRIBUTIONS--A HANDBOOK FOR C STUDENTS AND PRACTITIONERS, 1975, C PAGES 120-123. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING LABORATORY (205.03) C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE: 301-921-2315 C ORIGINAL VERSION--OCTOBER 1975. C UPDATED --NOVEMBER 1975. C C--------------------------------------------------------------------- C DOUBLE PRECISION PI DOUBLE PRECISION SQRT2 DOUBLE PRECISION DP DOUBLE PRECISION DNU DOUBLE PRECISION TERM1,TERM2,TERM3,TERM4,TERM5 DOUBLE PRECISION DPPFN DOUBLE PRECISION DPPF,DCON,DARG,Z,S,C DOUBLE PRECISION B21 DOUBLE PRECISION B31,B32,B33,B34 DOUBLE PRECISION B41,B42,B43,B44,B45 DOUBLE PRECISION B51,B52,B53,B54,B55,B56 DOUBLE PRECISION D1,D3,D5,D7,D9 DATA PI/3.14159265358979D0/ DATA SQRT2/1.414213562D0/ DATA B21/0.25D0/ DATA B31,B32,B33,B34/0.01041666666667D0,5.0D0,16.0D0,3.0D0/ DATA B41,B42,B43,B44,B45/0.00260416666667D0,3.0D0,19.0D0,17.0D0, 1 -15.0D0/ DATA B51,B52,B53,B54,B55,B56/0.00001085069444D0,79.0D0,776.0D0, 1 1482.0D0,-1920.0D0,-945.0D0/ C IPR=6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(P.LE.0.0.OR.P.GE.1.0)GOTO50 GOTO90 50 WRITE(IPR,1) WRITE(IPR,46)P RETURN 90 CONTINUE 1 FORMAT(1H ,115H***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO THE 1 TPPF SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****) 46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****) C C-----START POINT----------------------------------------------------- C DNU=NU DP=P MAXIT=5 C IF(NU.GE.3)GOTO250 IF(NU.EQ.1)GOTO100 IF(NU.EQ.2)GOTO200 WRITE(IPR,105) 105 FORMAT(1H ,33HINTERNAL ERROR IN TPPF SUBROUTINE) PPF=0.0 RETURN C C TREAT THE NU = 1 (CAUCHY) CASE C 100 DARG=PI*DP PPF=-DCOS(DARG)/DSIN(DARG) RETURN C C TREAT THE NU = 2 CASE C 200 TERM1=SQRT2/2.0D0 TERM2=2.0D0*DP-1.0D0 TERM3=DSQRT(DP*(1.0D0-DP)) PPF=TERM1*TERM2/TERM3 RETURN C C TREAT THE NU GREATER THAN OR EQUAL TO 3 CASE C 250 CALL NORPPF(P,PPFN) DPPFN=PPFN D1=DPPFN D3=DPPFN**3 D5=DPPFN**5 D7=DPPFN**7 D9=DPPFN**9 TERM1=D1 TERM2=B21*(D3+D1)/DNU TERM3=B31*(B32*D5+B33*D3+B34*D1)/(DNU**2) TERM4=B41*(B42*D7+B43*D5+B44*D3+B45*D1)/(DNU**3) TERM5=B51*(B52*D9+B53*D7+B54*D5+B55*D3+B56*D1)/(DNU**4) DPPF=TERM1+TERM2+TERM3+TERM4+TERM5 PPF=DPPF IF(NU.GE.7)RETURN IF(NU.EQ.3)GOTO300 IF(NU.EQ.4)GOTO400 IF(NU.EQ.5)GOTO500 IF(NU.EQ.6)GOTO600 RETURN C C AUGMENT THE RESULTS FOR THE NU = 3 CASE C 300 DCON=PI*(DP-0.5D0) DARG=DPPF/DSQRT(DNU) Z=DATAN(DARG) DO350IPASS=1,MAXIT S=DSIN(Z) C=DCOS(Z) Z=Z-(Z+S*C-DCON)/(2.0D0*C*C) 350 CONTINUE PPF=DSQRT(DNU)*S/C RETURN C C AUGMENT THE RESULTS FOR THE NU = 4 CASE C 400 DCON=2.0D0*(DP-0.5D0) DARG=DPPF/DSQRT(DNU) Z=DATAN(DARG) DO450IPASS=1,MAXIT S=DSIN(Z) C=DCOS(Z) Z=Z-((1.0D0+0.5D0*C*C)*S-DCON)/(1.5D0*C*C*C) 450 CONTINUE PPF=DSQRT(DNU)*S/C RETURN C C AUGMENT THE RESULTS FOR THE NU = 5 CASE C 500 DCON=PI*(DP-0.5D0) DARG=DPPF/DSQRT(DNU) Z=DATAN(DARG) DO550IPASS=1,MAXIT S=DSIN(Z) C=DCOS(Z) Z=Z-(Z+(C+(2.0D0/3.0D0)*C*C*C)*S-DCON)/((8.0D0/3.0D0)*C**4) 550 CONTINUE PPF=DSQRT(DNU)*S/C RETURN C C AUGMENT THE RESULTS FOR THE NU = 6 CASE C 600 DCON=2.0D0*(DP-0.5D0) DARG=DPPF/DSQRT(DNU) Z=DATAN(DARG) DO650IPASS=1,MAXIT S=DSIN(Z) C=DCOS(Z) Z=Z-((1.0D0+0.5D0*C*C+0.375D0*C**4)*S-DCON)/((15.0D0/8.0D0)*C**5) 650 CONTINUE PPF=DSQRT(DNU)*S/C RETURN C END