SUBROUTINE FCDF(X,NU1,NU2,CDF) C C PURPOSE--THIS SUBROUTINE COMPUTES THE CUMULATIVE DISTRIBUTION C FUNCTION VALUE FOR THE F DISTRIBUTION C WITH INTEGER DEGREES OF FREEDOM C PARAMETERS = NU1 AND NU2. C THIS DISTRIBUTION IS DEFINED FOR ALL NON-NEGATIVE X. C THE PROBABILITY DENSITY FUNCTION IS GIVEN C IN THE REFERENCES BELOW. C INPUT ARGUMENTS--X = THE SINGLE PRECISION VALUE AT C WHICH THE CUMULATIVE DISTRIBUTION C FUNCTION IS TO BE EVALUATED. C X SHOULD BE NON-NEGATIVE. C --NU1 = THE INTEGER DEGREES OF FREEDOM C FOR THE NUMERATOR OF THE F RATIO. C NU1 SHOULD BE POSITIVE. C --NU2 = THE INTEGER DEGREES OF FREEDOM C FOR THE DENOMINATOR OF THE F RATIO. C NU2 SHOULD BE POSITIVE. C OUTPUT ARGUMENTS--CDF = THE SINGLE PRECISION CUMULATIVE C DISTRIBUTION FUNCTION VALUE. C OUTPUT--THE SINGLE PRECISION CUMULATIVE DISTRIBUTION C FUNCTION VALUE CDF FOR THE F DISTRIBUTION C WITH DEGREES OF FREEDOM C PARAMETERS = NU1 AND NU2. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--X SHOULD BE NON-NEGATIVE. C --NU1 SHOULD BE A POSITIVE INTEGER VARIABLE. C --NU2 SHOULD BE A POSITIVE INTEGER VARIABLE. C OTHER DATAPAC SUBROUTINES NEEDED--NORCDF,CHSCDF. C FORTRAN LIBRARY SUBROUTINES NEEDED--DSQRT, DATAN. C MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION. C LANGUAGE--ANSI FORTRAN. C REFERENCES--NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS C SERIES 55, 1964, PAGES 946-947, C FORMULAE 26.6.4, 26.6.5, 26.6.8, AND 26.6.15. C --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE C DISTRIBUTIONS--2, 1970, PAGE 83, FORMULA 20, C AND PAGE 84, THIRD FORMULA. C --PAULSON, AN APPROXIMATE NORMAILIZATION C OF THE ANALYSIS OF VARIANCE DISTRIBUTION, C ANNALS OF MATHEMATICAL STATISTICS, 1942, C NUMBER 13, PAGES 233-135. C --SCHEFFE AND TUKEY, A FORMULA FOR SAMPLE SIZES C FOR POPULATION TOLERANCE LIMITS, 1944, C NUMBER 15, PAGE 217. 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--AUGUST 1972. C UPDATED --SEPTEMBER 1975. C UPDATED --NOVEMBER 1975. C UPDATED --OCTOBER 1976. C C--------------------------------------------------------------------- C DOUBLE PRECISION DX,PI,ANU1,ANU2,Z,SUM,TERM,AI,COEF1,COEF2,ARG DOUBLE PRECISION COEF DOUBLE PRECISION THETA,SINTH,COSTH,A,B DOUBLE PRECISION DSQRT,DATAN DOUBLE PRECISION DFACT1,DFACT2,DNUM,DDEN DOUBLE PRECISION DPOW1,DPOW2 DOUBLE PRECISION DNU1,DNU2 DOUBLE PRECISION TERM1,TERM2,TERM3 DATA PI/3.14159265358979D0/ DATA DPOW1,DPOW2/0.33333333333333D0,0.66666666666667D0/ DATA NUCUT1,NUCUT2/100,1000/ C IPR=6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(NU1.LE.0)GOTO50 IF(NU2.LE.0)GOTO55 IF(X.LT.0.0)GOTO60 GOTO90 50 WRITE(IPR,15) WRITE(IPR,47)NU1 CDF=0.0 RETURN 55 WRITE(IPR,23) WRITE(IPR,47)NU2 CDF=0.0 RETURN 60 WRITE(IPR,4) WRITE(IPR,46)X CDF=0.0 RETURN 90 CONTINUE 4 FORMAT(1H , 96H***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUME 1NT TO THE FCDF SUBROUTINE IS NEGATIVE *****) 15 FORMAT(1H , 91H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 FCDF SUBROUTINE IS NON-POSITIVE *****) 23 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD INPUT ARGUMENT TO THE 1 FCDF 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 DX=X M=NU1 N=NU2 ANU1=NU1 ANU2=NU2 DNU1=NU1 DNU2=NU2 C C IF X IS NON-POSITIVE, SET CDF = 0.0 AND RETURN. C IF NU2 IS 5 THROUGH 9 AND X IS MORE THAN 3000 C STANDARD DEVIATIONS BELOW THE MEAN, C SET CDF = 0.0 AND RETURN. C IF NU2 IS 10 OR LARGER AND X IS MORE THAN 150 C STANDARD DEVIATIONS BELOW THE MEAN, C SET CDF = 0.0 AND RETURN. C IF NU2 IS 5 THROUGH 9 AND X IS MORE THAN 3000 C STANDARD DEVIATIONS ABOVE THE MEAN, C SET CDF = 1.0 AND RETURN. C IF NU2 IS 10 OR LARGER AND X IS MORE THAN 150 C STANDARD DEVIATIONS ABOVE THE MEAN, C SET CDF = 1.0 AND RETURN. C IF(X.LE.0.0)GOTO105 IF(NU2.LE.4)GOTO109 T1=2.0/ANU1 T2=ANU2/(ANU2-2.0) T3=(ANU1+ANU2-2.0)/(ANU2-4.0) AMEAN=T2 SD=SQRT(T1*T2*T2*T3) ZRATIO=(X-AMEAN)/SD IF(NU2.LT.10.AND.ZRATIO.LT.-3000.0)GOTO105 IF(NU2.GE.10.AND.ZRATIO.LT.-150.0)GOTO105 IF(NU2.LT.10.AND.ZRATIO.GT.3000.0)GOTO107 IF(NU2.GE.10.AND.ZRATIO.GT.150.0)GOTO107 GOTO109 105 CDF=0.0 RETURN 107 CDF=1.0 RETURN 109 CONTINUE C C DISTINGUISH BETWEEN 6 SEPARATE REGIONS C OF THE (NU1,NU2) SPACE. C BRANCH TO THE PROPER COMPUTATIONAL METHOD C DEPENDING ON THE REGION. C NUCUT1 HAS THE VALUE 100. C NUCUT2 HAS THE VALUE 1000. C IF(NU1.LT.NUCUT2.AND.NU2.LT.NUCUT2)GOTO1000 IF(NU1.GE.NUCUT2.AND.NU2.GE.NUCUT2)GOTO2000 IF(NU1.LT.NUCUT1.AND.NU2.GE.NUCUT2)GOTO3000 IF(NU1.GE.NUCUT1.AND.NU2.GE.NUCUT2)GOTO2000 IF(NU1.GE.NUCUT2.AND.NU2.LT.NUCUT1)GOTO5000 IF(NU1.GE.NUCUT2.AND.NU2.GE.NUCUT1)GOTO2000 IBRAN=5 WRITE(IPR,99)IBRAN 99 FORMAT(1H ,42H*****INTERNAL ERROR IN FCDF SUBROUTINE--, 146HIMPOSSIBLE BRANCH CONDITION AT BRANCH POINT = ,I8) RETURN C C TREAT THE CASE WHEN NU1 AND NU2 C ARE BOTH SMALL OR MODERATE C (THAT IS, BOTH ARE SMALLER THAN 1000). C METHOD UTILIZED--EXACT FINITE SUM C (SEE AMS 55, PAGE 946, FORMULAE 26.6.4, 26.6.5, C AND 26.6.8). C 1000 CONTINUE Z=ANU2/(ANU2+ANU1*DX) IFLAG1=NU1-2*(NU1/2) IFLAG2=NU2-2*(NU2/2) IF(IFLAG1.EQ.0)GOTO120 IF(IFLAG2.EQ.0)GOTO150 GOTO250 C C DO THE NU1 EVEN AND NU2 EVEN OR ODD CASE C 120 SUM=0.0D0 TERM=1.0D0 IMAX=(M-2)/2 IF(IMAX.LE.0)GOTO110 DO100I=1,IMAX AI=I COEF1=2.0D0*(AI-1.0D0) COEF2=2.0D0*AI TERM=TERM*((ANU2+COEF1)/COEF2)*(1.0D0-Z) SUM=SUM+TERM 100 CONTINUE C 110 SUM=SUM+1.0D0 SUM=(Z**(ANU2/2.0D0))*SUM CDF=1.0D0-SUM RETURN C C DO THE NU1 ODD AND NU2 EVEN CASE C 150 SUM=0.0D0 TERM=1.0D0 IMAX=(N-2)/2 IF(IMAX.LE.0)GOTO210 DO200I=1,IMAX AI=I COEF1=2.0D0*(AI-1.0D0) COEF2=2.0D0*AI TERM=TERM*((ANU1+COEF1)/COEF2)*Z SUM=SUM+TERM 200 CONTINUE C 210 SUM=SUM+1.0D0 CDF=((1.0D0-Z)**(ANU1/2.0D0))*SUM RETURN C C DO THE NU1 ODD AND NU2 ODD CASE C 250 SUM=0.0D0 TERM=1.0D0 ARG=DSQRT((ANU1/ANU2)*DX) THETA=DATAN(ARG) SINTH=ARG/DSQRT(1.0D0+ARG*ARG) COSTH=1.0D0/DSQRT(1.0D0+ARG*ARG) IF(N.EQ.1)GOTO320 IF(N.EQ.3)GOTO310 IMAX=N-2 DO300I=3,IMAX,2 AI=I COEF1=AI-1.0D0 COEF2=AI TERM=TERM*(COEF1/COEF2)*(COSTH*COSTH) SUM=SUM+TERM 300 CONTINUE C 310 SUM=SUM+1.0D0 SUM=SUM*SINTH*COSTH C 320 A=(2.0D0/PI)*(THETA+SUM) 350 SUM=0.0D0 TERM=1.0D0 IF(M.EQ.1)B=0.0D0 IF(M.EQ.1)GOTO450 IF(M.EQ.3)GOTO410 IMAX=M-3 DO400I=1,IMAX,2 AI=I COEF1=AI COEF2=AI+2.0D0 TERM=TERM*((ANU2+COEF1)/COEF2)*(SINTH*SINTH) SUM=SUM+TERM 400 CONTINUE C 410 SUM=SUM+1.0D0 SUM=SUM*SINTH*(COSTH**N) COEF=1.0D0 IEVODD=N-2*(N/2) IMIN=3 IF(IEVODD.EQ.0)IMIN=2 IF(IMIN.GT.N)GOTO420 DO430I=IMIN,N,2 AI=I COEF=((AI-1.0D0)/AI)*COEF 430 CONTINUE C 420 COEF=COEF*ANU2 IF(IEVODD.EQ.0)GOTO440 COEF=COEF*(2.0D0/PI) C 440 B=COEF*SUM C 450 CDF=A-B RETURN C C TREAT THE CASE WHEN NU1 AND NU2 C ARE BOTH LARGE C (THAT IS, BOTH ARE EQUAL TO OR LARGER THAN 1000); C OR WHEN NU1 IS MODERATE AND NU2 IS LARGE C (THAT IS, WHEN NU1 IS EQUAL TO OR GREATER THAN 100 C BUT SMALLER THAN 1000, C AND NU2 IS EQUAL TO OR LARGER THAN 1000); C OR WHEN NU2 IS MODERATE AND NU1 IS LARGE C (THAT IS WHEN NU2 IS EQUAL TO OR GREATER THAN 100 C BUT SMALLER THAN 1000, C AND NU1 IS EQUAL TO OR LARGER THAN 1000). C METHOD UTILIZED--PAULSON APPROXIMATION C (SEE AMS 55, PAGE 947, FORMULA 26.6.15). C 2000 CONTINUE DFACT1=1.0D0/(4.5D0*DNU1) DFACT2=1.0D0/(4.5D0*DNU2) DNUM=((1.0D0-DFACT2)*(DX**DPOW1))-(1.0D0-DFACT1) DDEN=DSQRT((DFACT2*(DX**DPOW2))+DFACT1) U=DNUM/DDEN CALL NORCDF(U,GCDF) CDF=GCDF RETURN C C TREAT THE CASE WHEN NU1 IS SMALL C AND NU2 IS LARGE C (THAT IS, WHEN NU1 IS SMALLER THAN 100, C AND NU2 IS EQUAL TO OR LARGER THAN 1000). C METHOD UTILIZED--SHEFFE-TUKEY APPROXIMATION C (SEE JOHNSON AND KOTZ, VOLUME 2, PAGE 84, THIRD FORMULA). C 3000 CONTINUE TERM1=DNU1 TERM2=(DNU1/DNU2)*(0.5D0*DNU1-1.0D0) TERM3=-(DNU1/DNU2)*0.5D0 U=(TERM1+TERM2)/((1.0D0/DX)-TERM3) CALL CHSCDF(U,NU1,CCDF) CDF=CCDF RETURN C C TREAT THE CASE WHEN NU2 IS SMALL C AND NU1 IS LARGE C (THAT IS, WHEN NU2 IS SMALLER THAN 100, C AND NU1 IS EQUAL TO OR LARGER THAN 1000). C METHOD UTILIZED--SHEFFE-TUKEY APPROXIMATION C (SEE JOHNSON AND KOTZ, VOLUME 2, PAGE 84, THIRD FORMULA). C 5000 CONTINUE TERM1=DNU2 TERM2=(DNU2/DNU1)*(0.5D0*DNU2-1.0D0) TERM3=-(DNU2/DNU1)*0.5D0 U=(TERM1+TERM2)/(DX-TERM3) CALL CHSCDF(U,NU2,CCDF) CDF=1.0-CCDF RETURN C END