SUBROUTINE BINCDF(X,P,N,CDF) C C PURPOSE--THIS SUBROUTINE COMPUTES THE CUMULATIVE DISTRIBUTION C FUNCTION VALUE AT THE SINGLE PRECISION VALUE X C FOR THE BINOMIAL DISTRIBUTION C WITH SINGLE PRECISION 'BERNOULLI PROBABILITY' C PARAMETER = P, C AND INTEGER 'NUMBER OF BERNOULLI TRIALS' C PARAMETER = N. C THE BINOMIAL DISTRIBUTION USED C HEREIN HAS MEAN = N*P C AND STANDARD DEVIATION = SQRT(N*P*(1-P)). C THIS DISTRIBUTION IS DEFINED FOR ALL C DISCRETE INTEGER X BETWEEN 0 (INCLUSIVELY) C AND N (INCLUSIVELY). C THIS DISTRIBUTION HAS THE PROBABILITY FUNCTION C F(X) = C(N,X) * P**X * (1-P)**(N-X). C WHERE C(N,X) IS THE COMBINATORIAL FUNCTION C EQUALING THE NUMBER OF COMBINATIONS OF N ITEMS C TAKEN X AT A TIME. C THE BINOMIAL DISTRIBUTION IS THE C DISTRIBUTION OF THE NUMBER OF C SUCCESSES IN N BERNOULLI (0,1) C TRIALS WHERE THE PROBABILITY OF SUCCESS C IN A SINGLE TRIAL = P. C INPUT ARGUMENTS--X = THE SINGLE PRECISION VALUE C AT WHICH THE CUMULATIVE DISTRIBUTION C FUNCTION IS TO BE EVALUATED. C X SHOULD BE INTEGRAL-VALUED, C AND BETWEEN 0.0 (INCLUSIVELY) C AND N (INCLUSIVELY). C --P = THE SINGLE PRECISION VALUE C OF THE 'BERNOULLI PROBABILITY' C PARAMETER FOR THE BINOMIAL C DISTRIBUTION. C P SHOULD BE BETWEEN C 0.0 (EXCLUSIVELY) AND C 1.0 (EXCLUSIVELY). C --N = THE INTEGER VALUE C OF THE 'NUMBER OF BERNOULLI TRIALS' C PARAMETER. C N SHOULD BE A POSITIVE INTEGER. C OUTPUT ARGUMENTS--CDF = THE SINGLE PRECISION CUMULATIVE C DISTRIBUTION FUNCTION VALUE. C OUTPUT--THE SINGLE PRECISION CUMULATIVE DISTRIBUTION C FUNCTION VALUE CDF C FOR THE BINOMIAL DISTRIBUTION C WITH 'BERNOULLI PROBABILITY' PARAMETER = P C AND 'NUMBER OF BERNOULLI TRIALS' PARAMETER = N. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--X SHOULD BE INTEGRAL-VALUED, C AND BETWEEN 0.0 (INCLUSIVELY) C AND N (INCLUSIVELY). C --P SHOULD BE BETWEEN 0.0 (EXCLUSIVELY) C AND 1.0 (EXCLUSIVELY). C --N SHOULD BE A POSITIVE INTEGER. C OTHER DATAPAC SUBROUTINES NEEDED--NONE. C FORTRAN LIBRARY SUBROUTINES NEEDED--DSQRT, DATAN. C MODE OF INTERNAL OPERATIONS--DOUBLE PRECISION. C LANGUAGE--ANSI FORTRAN. C COMMENT--NOTE THAT EVEN THOUGH THE INPUT C TO THIS CUMULATIVE C DISTRIBUTION FUNCTION SUBROUTINE C FOR THIS DISCRETE DISTRIBUTION C SHOULD (UNDER NORMAL CIRCUMSTANCES) BE A C DISCRETE INTEGER VALUE, C THE INPUT VARIABLE X IS SINGLE C PRECISION IN MODE. C X HAS BEEN SPECIFIED AS SINGLE C PRECISION SO AS TO CONFORM WITH THE DATAPAC C CONVENTION THAT ALL INPUT ****DATA**** C (AS OPPOSED TO SAMPLE SIZE, FOR EXAMPLE) C VARIABLES TO ALL C DATAPAC SUBROUTINES ARE SINGLE PRECISION. C THIS CONVENTION IS BASED ON THE BELIEF THAT C 1) A MIXTURE OF MODES (FLOATING POINT C VERSUS INTEGER) IS INCONSISTENT AND C AN UNNECESSARY COMPLICATION C IN A DATA ANALYSIS; AND C 2) FLOATING POINT MACHINE ARITHMETIC C (AS OPPOSED TO INTEGER ARITHMETIC) C IS THE MORE NATURAL MODE FOR DOING C DATA ANALYSIS. C REFERENCES--HASTINGS AND PEACOCK, STATISTICAL C DISTRIBUTIONS--A HANDBOOK FOR C STUDENTS AND PRACTITIONERS, 1975, C PAGE 38. C --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS C SERIES 55, 1964, PAGE 945, FORMULAE 26.5.24 AND C 26.5.28, AND PAGE 929. C --JOHNSON AND KOTZ, DISCRETE C DISTRIBUTIONS, 1969, PAGES 50-86, C ESPECIALLY PAGES 63-64. C --FELLER, AN INTRODUCTION TO PROBABILITY C THEORY AND ITS APPLICATIONS, VOLUME 1, C EDITION 2, 1957, PAGES 135-142. C --KENDALL AND STUART, THE ADVANCED THEORY OF C STATISTICS, VOLUME 1, EDITION 2, 1963, PAGES 120-125. C --MOOD AND GRABLE, INTRODUCTION TO THE THEORY C OF STATISTICS, EDITION 2, 1963, PAGES 64-69. C --OWEN, HANDBOOK OF STATISTICAL C TABLES, 1962, PAGES 264-272. 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--NOVEMBER 1975. C UPDATED --MAY 1977. 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 DATA PI/3.14159265358979D0/ C IPR=6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C AN=N IF(P.LT.0.0.OR.P.GT.1.0)GOTO50 IF(N.LT.1)GOTO55 IF(X.LT.0.0.OR.X.GT.AN)GOTO60 INTX=X+0.0001 FINTX=INTX DEL=X-FINTX IF(DEL.LT.0.0)DEL=-DEL IF(DEL.GT.0.001)GOTO65 GOTO90 50 WRITE(IPR,11) WRITE(IPR,46)P CDF=0.0 RETURN 55 WRITE(IPR,25) WRITE(IPR,47)N CDF=0.0 RETURN 60 WRITE(IPR,4)N WRITE(IPR,46)X IF(X.LT.0.0)CDF=0.0 IF(X.GT.AN)CDF=1.0 RETURN 65 WRITE(IPR,5) WRITE(IPR,46)X 90 CONTINUE 4 FORMAT(1H ,111H***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUME 1NT TO THE BINCDF SUBROUTINE IS OUTSIDE THE USUAL (0,N) = (0,,I7, 1 11H,INTERVAL *) 5 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUME 1NT TO THE BINCDF SUBROUTINE IS NON-INTEGRAL *****) 11 FORMAT(1H ,115H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 BINCDF SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****) 25 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD INPUT ARGUMENT TO THE 1 BINCDF 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 TREAT IMMEDIATELY THE SPECIAL CASE OF X = N, C IN WHICH CASE CDF = 1.0. C ALSO TREAT IMMEDIATELY THE SPECIAL CASE OF P = 0.0 C IN WHICH CASE CDF = 1.0 FOR ALL X. C THIRDLY, TREAT THE SPECIAL CASE IN WHICH P = 1.0 C IN WHICH CASE CDF = 0.0 FOR ALL X SMALLER THAN N C AND CDF = 1.0 FOR ALL X EQUAL TO OR LARGER C THAN N. C INTX=X+0.0001 CDF=1.0 IF(INTX.EQ.N)RETURN IF(P.EQ.0.0)RETURN IF(P.EQ.1.0.AND.INTX.GE.N)RETURN IF(P.EQ.1.0.AND.INTX.LT.N)CDF=0.0 IF(P.EQ.1.0.AND.INTX.LT.N)RETURN C C EXPRESS THE BINOMIAL CUMULATIVE DISTRIBUTION C FUNCTION IN TERMS OF THE EQUIVALENT F C CUMULATIVE DISTRIBUTION FUNCTION, C AND THEN EVALUATE THE LATTER. C AN=N DX=(P/(1.0-P))*((AN-X)/(X+1.0)) NU1=2.0*(X+1.0)+0.1 NU2=2.0*(AN-X)+0.1 ANU1=NU1 ANU2=NU2 Z=ANU2/(ANU2+ANU1*DX) C C DETERMINE IF NU1 AND NU2 ARE EVEN OR ODD C 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=(NU1-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=SUM RETURN C C DO THE NU1 ODD AND NU2 EVEN CASE C 150 SUM=0.0D0 TERM=1.0D0 IMAX=(NU2-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-((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(NU2.EQ.1)GOTO320 IF(NU2.EQ.3)GOTO310 IMAX=NU2-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) SUM=0.0D0 TERM=1.0D0 IF(NU1.EQ.1)B=0.0D0 IF(NU1.EQ.1)GOTO450 IF(NU1.EQ.3)GOTO410 IMAX=NU1-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=NU2-2*(NU2/2) IMIN=3 IF(IEVODD.EQ.0)IMIN=2 IF(IMIN.GT.NU2)GOTO420 DO430I=IMIN,NU2,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=1.0D0-(A-B) RETURN C END