SUBROUTINE NBCDF(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 NEGATIVE BINOMIAL DISTRIBUTION C WITH SINGLE PRECISION 'BERNOULLI PROBABILITY' C PARAMETER = P, C AND INTEGER 'NUMBER OF SUCCESSES IN BERNOULLI TRIALS' C PARAMETER = N. C THE NEGATIVE BINOMIAL DISTRIBUTION USED C HEREIN HAS MEAN = N*(1-P)/P C AND STANDARD DEVIATION = SQRT(N*(1-P)/(P*P))). C THIS DISTRIBUTION IS DEFINED FOR C ALL NON-NEGATIVE INTEGER X--X = 0, 1, 2, ... . C THIS DISTRIBUTION HAS THE PROBABILITY FUNCTION C F(X) = C(N+X-1,N) * P**N * (1-P)**X. C WHERE C(N+X-1,N) IS THE COMBINATORIAL FUNCTION C EQUALING THE NUMBER OF COMBINATIONS OF N+X-1 ITEMS C TAKEN N AT A TIME. C THE NEGATIVE BINOMIAL DISTRIBUTION IS THE C DISTRIBUTION OF THE NUMBER OF FAILURES C BEFORE OBTAINING N SUCCESSES IN AN C INDEFINITE SEQUENCE OF 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 NON-NEGATIVE AND C INTEGRAL-VALUED. C --P = THE SINGLE PRECISION VALUE C OF THE 'BERNOULLI PROBABILITY' C PARAMETER FOR THE NEGATIVE 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 SUCCESSES C IN BERNOULLI TRIALS' 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 NEGATIVE BINOMIAL DISTRIBUTION C WITH 'BERNOULLI PROBABILITY' PARAMETER = P C AND 'NUMBER OF SUCCESSES IN BERNOULLI TRIALS' C PARAMETER = N. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--X SHOULD BE NON-NEGATIVE AND INTEGRAL-VALUED. 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--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 122-142, C ESPECIALLY PAGE 127. C --HASTINGS AND PEACOCK, STATISTICAL C DISTRIBUTIONS--A HANDBOOK FOR C STUDENTS AND PRACTITIONERS, 1975, C PAGES 92-95. C --FELLER, AN INTRODUCTION TO PROBABILITY C THEORY AND ITS APPLICATIONS, VOLUME 1, C EDITION 2, 1957, PAGES 155-157, 210. C --KENDALL AND STUART, THE ADVANCED THEORY OF C STATISTICS, VOLUME 1, EDITION 2, 1963, PAGES 130-131. C --WILLIAMSON AND BRETHERTON, TABLES OF C THE NEGATIVE BINOMIAL PROBABILITY C DISTRIBUTION, 1963. C --OWEN, HANDBOOK OF STATISTICAL C TABLES, 1962, PAGE 304. 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 C--------------------------------------------------------------------- C DOUBLE PRECISION DX2,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.LE.0.0.OR.P.GE.1.0)GOTO50 IF(N.LT.1)GOTO55 IF(X.LT.0.0)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) WRITE(IPR,46)X IF(X.LT.0.0)CDF=0.0 RETURN 65 WRITE(IPR,5) WRITE(IPR,46)X 90 CONTINUE 4 FORMAT(1H , 96H***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUME 1NT TO THE NBCDF SUBROUTINE IS NEGATIVE *****) 5 FORMAT(1H ,100H***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUME 1NT TO THE NBCDF SUBROUTINE IS NON-INTEGRAL *****) 11 FORMAT(1H ,115H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 NBCDF SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****) 25 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD INPUT ARGUMENT TO THE 1 NBCDF 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 NEGATIVE BINOMIAL CUMULATIVE DISTRIBUTION C FUNCTION IN TERMS OF THE EQUIVALENT BINOMIAL C CUMULATIVE DISTRIBUTION FUNCTION, C AND THEN OPERATE ON THE LATTER. C INTX=X+0.0001 K=N-1 N2=N+INTX 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 AK=K AN2=N2 DX2=(P/(1.0-P))*((AN2-AK)/(AK+1.0)) NU1=2*(K+1) NU2=2*(N2-K) ANU1=NU1 ANU2=NU2 Z=ANU2/(ANU2+ANU1*DX2) 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=1.0D0-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-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)*DX2) 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) 350 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=A-B RETURN C END