SUBROUTINE NBPPF(P,PPAR,N,PPF) C C PURPOSE--THIS SUBROUTINE COMPUTES THE PERCENT POINT C FUNCTION VALUE AT THE SINGLE PRECISION VALUE P C FOR THE NEGATIVE BINOMIAL DISTRIBUTION C WITH SINGLE PRECISION 'BERNOULLI PROBABILITY' C PARAMETER = PPAR, C AND INTEGER 'NUMBER OF SUCCESSES IN BERNOULLI TRIALS' C PARAMETER = N. C THE NEGATIVE BINOMIAL DISTRIBUTION USED C HEREIN HAS MEAN = N*(1-PPAR)/PPAR C AND STANDARD DEVIATION = SQRT(N*(1-PPAR)/(PPAR*PPAR))). 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) * PPAR**N * (1-PPAR)**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 = PPAR. 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 --PPAR = THE SINGLE PRECISION VALUE C OF THE 'BERNOULLI PROBABILITY' C PARAMETER FOR THE NEGATIVE BINOMIAL C DISTRIBUTION. C PPAR 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--PPF = THE SINGLE PRECISION PERCENT C POINT FUNCTION VALUE. C OUTPUT--THE SINGLE PRECISION PERCENT POINT . C FUNCTION VALUE PPF C FOR THE NEGATIVE BINOMIAL DISTRIBUTION C WITH 'BERNOULLI PROBABILITY' PARAMETER = PPAR C AND 'NUMBER OF SUCCESSES IN BERNOULLI TRIALS' C PARAMETER = N. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--PPAR SHOULD BE BETWEEN 0.0 (EXCLUSIVELY) C AND 1.0 (EXCLUSIVELY). C --N SHOULD BE A POSITIVE INTEGER. C --P SHOULD BE BETWEEN 0.0 (INCLUSIVELY) C AND 1.0 (EXCLUSIVELY). C OTHER DATAPAC SUBROUTINES NEEDED--NORPPF, NBCDF. C FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, EXP, ALOG. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION AND DOUBLE PRECISION. C LANGUAGE--ANSI FORTRAN. C COMMENT--NOTE THAT EVEN THOUGH THE OUTPUT C FROM THIS DISCRETE DISTRIBUTION C PERCENT POINT FUNCTION C SUBROUTINE MUST NECESSARILY BE A C DISCRETE INTEGER VALUE, C THE OUTPUT VARIABLE PPF IS SINGLE C PRECISION IN MODE. C PPF HAS BEEN SPECIFIED AS SINGLE C PRECISION SO AS TO CONFORM WITH THE DATAPAC C CONVENTION THAT ALL OUTPUT VARIABLES FROM 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--JOHNSON AND KOTZ, DISCRETE C DISTRIBUTIONS, 1969, PAGES 122-142, C ESPECIALLY PAGE 127, FORMULA 22. C --HASTINGS AND PEACOCK, STATISTICAL C DISTRIBUTIONS--A HANDBOOK FOR C STUDENTS AND PRACTITIONERS, 1975, C PAGES 92-95. C --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS C SERIES 55, 1964, PAGE 929. 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 DPPAR C IPR=6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(P.LT.0.0.OR.P.GE.1.0)GOTO50 IF(PPAR.LE.0.0.OR.PPAR.GE.1.0)GOTO55 IF(N.LT.1)GOTO60 GOTO90 50 WRITE(IPR,1) WRITE(IPR,46)P PPF=0.0 RETURN 55 WRITE(IPR,11) WRITE(IPR,46)PPAR PPF=0.0 RETURN 60 WRITE(IPR,25) WRITE(IPR,47)N PPF=0.0 RETURN 90 CONTINUE 1 FORMAT(1H ,115H***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO THE 1 NBPPF SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****) 11 FORMAT(1H ,115H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 NBPPF SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****) 25 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD INPUT ARGUMENT TO THE 1 NBPPF 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 AN=N DPPAR=PPAR PPF=0.0 IX0=0 IX1=0 IX2=0 P0=0.0 P1=0.0 P2=0.0 C C TREAT CERTAIN SPECIAL CASES IMMEDIATELY-- C 1) P = 0.0 C 2) P = 0.5 AND PPAR = 0.5 C 3) PPF = 0 C IF(P.EQ.0.0)GOTO110 IF(P.EQ.0.5.AND.PPAR.EQ.0.5)GOTO130 PF0=DPPAR**N IF(P.LE.PF0)GOTO110 GOTO190 110 PPF=0.0 RETURN 130 PPF=N-1 RETURN 190 CONTINUE C C DETERMINE AN INITIAL APPROXIMATION TO THE NEGATIVE BINOMIAL C PERCENT POINT BY USE OF THE HYPERBOLIC ARCSIN C TRANSFORMATION OF THE NEGATIVE BINOMIAL C TO APPROXIMATE NORMALITY. C (SEE JOHNSON AND KOTZ, DISCRETE DISTRIBUTIONS, C PAGE 127, FORMULA 22). C AMEAN=AN*(1.0-PPAR)/PPAR SD=SQRT(AN*(1.0-PPAR)/(PPAR*PPAR)) ARG=SQRT((AMEAN+0.375)/(AN-0.75)) ARCSH=ALOG(ARG+SQRT(ARG*ARG+1.0)) YMEAN=(SQRT(AN-0.5))*ARCSH YSD=0.5 CALL NORPPF(P,ZPPF) YPPF=YMEAN+ZPPF*YSD ARG=YPPF/SQRT(AN-0.5) E=EXP(ARG) SINH=(E-1.0/E)/2.0 X2=-0.375+(AN-0.75)*SINH*SINH X2=X2+0.5 IX2=X2 C C CHECK AND MODIFY (IF NECESSARY) THIS INITIAL C ESTIMATE OF THE PERCENT POINT C TO ASSURE THAT IT BE NON-NEGATIVE. C IF(IX2.LT.0)IX2=0 C C DETERMINE UPPER AND LOWER BOUNDS ON THE DESIRED C PERCENT POINT BY ITERATING OUT (BOTH BELOW AND ABOVE) C FROM THE ORIGINAL APPROXIMATION AT STEPS C OF 1 STANDARD DEVIATION. C THE RESULTING BOUNDS WILL BE AT MOST C 1 STANDARD DEVIATION APART. C IX0=0 IX1=10**10 ISD=SD+1.0 X2=IX2 CALL NBCDF(X2,PPAR,N,P2) C IF(P2.LT.P)GOTO210 GOTO250 C 210 IX0=IX2 DO220I=1,100000 IX2=IX0+ISD IF(IX2.GE.IX1)GOTO275 X2=IX2 CALL NBCDF(X2,PPAR,N,P2) IF(P2.GE.P)GOTO230 IX0=IX2 220 CONTINUE WRITE(IPR,249) WRITE(IPR,222) GOTO950 230 IX1=IX2 GOTO275 C 250 IX1=IX2 DO260I=1,100000 IX2=IX1-ISD IF(IX2.LE.IX0)GOTO275 X2=IX2 CALL NBCDF(X2,PPAR,N,P2) IF(P2.LT.P)GOTO270 IX1=IX2 260 CONTINUE WRITE(IPR,249) WRITE(IPR,262) GOTO950 270 IX0=IX2 C 275 IF(IX0.EQ.IX1)GOTO280 GOTO295 280 IF(IX0.EQ.0)GOTO285 IF(IX0.EQ.N)GOTO290 WRITE(IPR,249) WRITE(IPR,282) GOTO950 285 IX1=IX1+1 GOTO295 290 IX0=IX0-1 295 CONTINUE C C COMPUTE NEGATIVE BINOMIAL PROBABILITIES FOR THE C DERIVED LOWER AND UPPER BOUNDS. C X0=IX0 X1=IX1 CALL NBCDF(X0,PPAR,N,P0) CALL NBCDF(X1,PPAR,N,P1) C C CHECK THE PROBABILITIES FOR PROPER ORDERING C IF(P0.LT.P.AND.P.LE.P1)GOTO490 IF(P0.EQ.P)GOTO410 IF(P1.EQ.P)GOTO420 IF(P0.GT.P1)GOTO430 IF(P0.GT.P)GOTO440 IF(P1.LT.P)GOTO450 WRITE(IPR,249) WRITE(IPR,401) GOTO950 410 PPF=IX0 RETURN 420 PPF=IX1 RETURN 430 WRITE(IPR,249) WRITE(IPR,431) GOTO950 440 WRITE(IPR,249) WRITE(IPR,441) GOTO950 450 WRITE(IPR,249) WRITE(IPR,451) GOTO950 490 CONTINUE C C THE STOPPING CRITERION IS THAT THE LOWER BOUND C AND UPPER BOUND ARE EXACTLY 1 UNIT APART. C CHECK TO SEE IF IX1 = IX0 + 1; C IF SO, THE ITERATIONS ARE COMPLETE; C IF NOT, THEN BISECT, COMPUTE PROBABILIIES, C CHECK PROBABILITIES, AND CONTINUE ITERATING C UNTIL IX1 = IX0 + 1. C 300 IX0P1=IX0+1 IF(IX1.EQ.IX0P1)GOTO690 IX2=(IX0+IX1)/2 IF(IX2.EQ.IX0)GOTO610 IF(IX2.EQ.IX1)GOTO620 X2=IX2 CALL NBCDF(X2,PPAR,N,P2) IF(P0.LT.P2.AND.P2.LT.P1)GOTO630 IF(P2.LE.P0)GOTO640 IF(P2.GE.P1)GOTO650 610 WRITE(IPR,249) WRITE(IPR,611) GOTO950 620 WRITE(IPR,249) WRITE(IPR,611) GOTO950 630 IF(P2.LE.P)GOTO635 IX1=IX2 P1=P2 GOTO300 635 IX0=IX2 P0=P2 GOTO300 640 WRITE(IPR,249) WRITE(IPR,641) GOTO950 650 WRITE(IPR,249) WRITE(IPR,651) GOTO950 690 PPF=IX1 IF(P0.EQ.P)PPF=IX0 RETURN C 950 WRITE(IPR,240)IX0,P0 WRITE(IPR,241)IX1,P1 WRITE(IPR,242)IX2,P2 WRITE(IPR,244)P WRITE(IPR,245)PPAR,N RETURN C 222 FORMAT(1H ,43HNO UPPER BOUND FOUND AFTER 10**7 ITERATIONS) 240 FORMAT(1H ,7HIX0 = ,I8,10X,5HP0 = ,F14.7) 241 FORMAT(1H ,7HIX1 = ,I8,10X,5HP1 = ,F14.7) 242 FORMAT(1H ,7HIX2 = ,I8,10X,5HP2 = ,F14.7) 244 FORMAT(1H ,7HP = ,F14.7) 245 FORMAT(1H ,7HPPAR = ,F14.7,10X,5HN = ,I8) 249 FORMAT(1H ,47H***** INTERNAL ERROR IN NBPPF SUBROUTINE *****) 262 FORMAT(1H ,43HNO LOWER BOUND FOUND AFTER 10**7 ITERATIONS) 282 FORMAT(1H ,31HLOWER AND UPPER BOUND IDENTICAL) 401 FORMAT(1H ,39HIMPOSSIBLE BRANCH CONDITION ENCOUNTERED) 431 FORMAT(1H ,42HLOWER BOUND PROBABILITY (P0) GREATER THAN , 1 28HUPPER BOUND PROBABILITY (P1)) 441 FORMAT(1H ,42HLOWER BOUND PROBABILITY (P0) GREATER THAN , 1 21HINPUT PROBABILITY (P)) 451 FORMAT(1H ,42HUPPER BOUND PROBABILITY (P1) LESS THAN , 1 21HINPUT PROBABILITY (P)) 611 FORMAT(1H ,39HBISECTION VALUE (X2) = LOWER BOUND (X0)) 621 FORMAT(1H ,39HBISECTION VALUE (X2) = UPPER BOUND (X1)) 641 FORMAT(1H ,33HBISECTION VALUE PROBABILITY (P2) , 1 38HLESS THAN LOWER BOUND PROBABILITY (P0)) 651 FORMAT(1H ,33HBISECTION VALUE PROBABILITY (P2) , 1 41HGREATER THAN UPPER BOUND PROBABILITY (P1)) C END