SUBROUTINE CHSPPF(P,NU,PPF) C C PURPOSE--THIS SUBROUTINE COMPUTES THE PERCENT POINT C FUNCTION VALUE FOR THE CHI-SQUARED DISTRIBUTION C WITH INTEGER DEGREES OF FREEDOM PARAMETER = NU. C THE CHI-SQUARED DISTRIBUTION USED C HEREIN IS DEFINED FOR ALL NON-NEGATIVE X, C AND ITS PROBABILITY DENSITY FUNCTION IS GIVEN C IN REFERENCES 2, 3, AND 4 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 (INCLUSIVELY) 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 CHI-SQUARED 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 (INCLUSIVELY) C AND 1.0 (EXCLUSIVELY). C OTHER DATAPAC SUBROUTINES NEEDED--NONE. C FORTRAN LIBRARY SUBROUTINES NEEDED--DEXP, DLOG. C MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION. C LANGUAGE--ANSI FORTRAN. C ACCURACY--(ON THE UNIVAC 1108, EXEC 8 SYSTEM AT NBS) C COMPARED TO THE KNOWN NU = 2 (EXPONENTIAL) C RESULTS, AGREEMENT WAS HAD OUT TO 6 SIGNIFICANT C DIGITS FOR ALL TESTED P IN THE RANGE P = .001 TO C P = .999. FOR P = .95 AND SMALLER, THE AGREEMENT C WAS EVEN BETTER--7 SIGNIFICANT DIGITS. C (NOTE THAT THE TABULATED VALUES GIVEN IN THE WILK, C GNANADESIKAN, AND HUYETT REFERENCE BELOW, PAGE 20, C ARE IN ERROR FOR AT LEAST THE GAMMA = 1 CASE-- C THE WORST DETECTED ERROR WAS AGREEMENT TO ONLY 3 C SIGNIFICANT DIGITS (IN THEIR 8 SIGNIFICANT DIGIT TABLE) C FOR P = .999.) C REFERENCES--WILK, GNANADESIKAN, AND HUYETT, 'PROBABILITY C PLOTS FOR THE GAMMA DISTRIBUTION', C TECHNOMETRICS, 1962, PAGES 1-15, C ESPECIALLY PAGES 3-5. C --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS C SERIES 55, 1964, PAGE 257, FORMULA 6.1.41, C AND PAGES 940-943. C --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE C DISTRIBUTIONS--1, 1970, PAGES 166-206. C --HASTINGS AND PEACOCK, STATISTICAL C DISTRIBUTIONS--A HANDBOOK FOR C STUDENTS AND PRACTITIONERS, 1975, C PAGES 46-51. 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--SEPTEMBER 1975. C UPDATED --NOVEMBER 1975. C C--------------------------------------------------------------------- C DOUBLE PRECISION DP,DGAMMA DOUBLE PRECISION Z,Z2,Z3,Z4,Z5,DEN,A,B,C,D,G DOUBLE PRECISION XMIN0,XMIN,AI,XMAX,DX,PCALC,XMID DOUBLE PRECISION XLOWER,XUPPER,XDEL DOUBLE PRECISION SUM,TERM,CUT1,CUT2,AJ,CUTOFF,T DOUBLE PRECISION DEXP,DLOG DIMENSION D(10) DATA C/ .918938533204672741D0/ DATA D(1),D(2),D(3),D(4),D(5) 1 /+.833333333333333333D-1,-.277777777777777778D-2, 1+.793650793650793651D-3,-.595238095238095238D-3,+.8417508417508417 151D-3/ DATA D(6),D(7),D(8),D(9),D(10) 1 /-.191752691752691753D-2,+.641025641025641025D-2,-.2955065359 147712418D-1,+.179644372368830573D0,-.139243221690590111D1/ C IPR=6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(P.LT.0.0.OR.P.GE.1.0)GOTO50 IF(NU.LT.1)GOTO55 GOTO90 50 WRITE(IPR,1) WRITE(IPR,46)P PPF=0.0 RETURN 55 WRITE(IPR,15) WRITE(IPR,47)NU PPF=0.0 RETURN 90 CONTINUE 1 FORMAT(1H ,115H***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO THE 1 CHSPPF SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****) 15 FORMAT(1H , 91H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 CHSPPF SUBROUTINE IS NON-POSITIVE *****) 46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****) 47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8 ,6H *****) C C-----START POINT----------------------------------------------------- C C EXPRESS THE CHI-SQUARED DISTRIBUTION PERCENT POINT C FUNCTION IN TERMS OF THE EQUIVALENT GAMMA C DISTRIBUTION PERCENT POINT FUNCTION, C AND THEN EVALUATE THE LATTER. C ANU=NU GAMMA=ANU/2.0 DP=P DNU=NU DGAMMA=DNU/2.0D0 MAXIT=10000 C C COMPUTE THE GAMMA FUNCTION USING THE ALGORITHM IN THE C NBS APPLIED MATHEMATICS SERIES REFERENCE. C THIS GAMMA FUNCTION NEED BE CALCULATED ONLY ONCE. C IT IS USED IN THE CALCULATION OF THE CDF BASED ON C THE TENTATIVE VALUE OF THE PPF IN THE ITERATION. C Z=DGAMMA DEN=1.0D0 150 IF(Z.GE.10.0D0)GOTO160 DEN=DEN*Z Z=Z+1.0D0 GOTO150 160 Z2=Z*Z Z3=Z*Z2 Z4=Z2*Z2 Z5=Z2*Z3 A=(Z-0.5D0)*DLOG(Z)-Z+C B=D(1)/Z+D(2)/Z3+D(3)/Z5+D(4)/(Z2*Z5)+D(5)/(Z4*Z5)+ 1D(6)/(Z*Z5*Z5)+D(7)/(Z3*Z5*Z5)+D(8)/(Z5*Z5*Z5)+D(9)/(Z2*Z5*Z5*Z5) G=DEXP(A+B)/DEN C C DETERMINE LOWER AND UPPER LIMITS ON THE DESIRED 100P C PERCENT POINT. C ILOOP=1 XMIN0=(DP*DGAMMA*G)**(1.0D0/DGAMMA) XMIN=XMIN0 ICOUNT=1 350 AI=ICOUNT XMAX=AI*XMIN0 DX=XMAX GOTO1000 360 IF(PCALC.GE.DP)GOTO370 XMIN=XMAX ICOUNT=ICOUNT+1 IF(ICOUNT.LE.30000)GOTO350 370 XMID=(XMIN+XMAX)/2.0D0 C C NOW ITERATE BY BISECTION UNTIL THE DESIRED ACCURACY IS ACHIEVED. C ILOOP=2 XLOWER=XMIN XUPPER=XMAX ICOUNT=0 550 DX=XMID GOTO1000 560 IF(PCALC.EQ.DP)GOTO570 IF(PCALC.GT.DP)GOTO580 XLOWER=XMID XMID=(XMID+XUPPER)/2.0D0 GOTO590 580 XUPPER=XMID XMID=(XMID+XLOWER)/2.0D0 590 XDEL=XMID-XLOWER IF(XDEL.LT.0.0D0)XDEL=-XDEL ICOUNT=ICOUNT+1 IF(XDEL.LT.0.0000000001D0.OR.ICOUNT.GT.100)GOTO570 GOTO550 570 PPF=2.0D0*XMID RETURN C C******************************************************************** C THIS SECTION BELOW IS LOGICALLY SEPARATE FROM THE ABOVE. C THIS SECTION COMPUTES A CDF VALUE FOR ANY GIVEN TENTATIVE C PERCENT POINT X VALUE AS DEFINED IN EITHER OF THE 2 C ITERATION LOOPS IN THE ABOVE CODE. C C COMPUTE T-SUB-Q AS DEFINED ON PAGE 4 OF THE WILK, GNANADESIKAN, C AND HUYETT REFERENCE C 1000 SUM=1.0D0/DGAMMA TERM=1.0D0/DGAMMA CUT1=DX-DGAMMA CUT2=DX*10000000000.0D0 DO700J=1,MAXIT AJ=J TERM=DX*TERM/(DGAMMA+AJ) SUM=SUM+TERM CUTOFF=CUT1+(CUT2*TERM/SUM) IF(AJ.GT.CUTOFF)GOTO750 700 CONTINUE WRITE(IPR,705)MAXIT WRITE(IPR,706)P WRITE(IPR,707)NU WRITE(IPR,708) PPF=0.0 RETURN C 750 T=SUM PCALC=(DX**DGAMMA)*(DEXP(-DX))*T/G IF(ILOOP.EQ.1)GOTO360 GOTO560 C 705 FORMAT(1H ,48H*****ERROR IN INTERNAL OPERATIONS IN THE CHSPPF , 1 45HSUBROUTINE--THE NUMBER OF ITERATIONS EXCEEDS ,I7) 706 FORMAT(1H ,33H THE INPUT VALUE OF P IS ,E15.8) 707 FORMAT(1H ,33H THE INPUT VALUE OF NU IS ,I8) 708 FORMAT(1H ,48H THE OUTPUT VALUE OF PPF HAS BEEN SET TO 0.0) C END