SUBROUTINE POIPPF(P,ALAMBA,PPF) C C PURPOSE--THIS SUBROUTINE COMPUTES THE PERCENT POINT C FUNCTION VALUE AT THE SINGLE PRECISION VALUE P C FOR THE POISSON DISTRIBUTION C WITH SINGLE PRECISION C TAIL LENGTH PARAMETER = ALAMBA. C THE POISSON DISTRIBUTION USED C HEREIN HAS MEAN = ALAMBA C AND STANDARD DEVIATION = SQRT(ALAMBA). C THIS DISTRIBUTION IS DEFINED FOR C ALL DISCRETE NON-NEGATIVE INTEGER X--X = 0, 1, 2, ... . C THIS DISTRIBUTION HAS THE PROBABILITY FUNCTION C F(X) = EXP(-ALAMBA) * ALAMBA**X / X!. C THE POISSON DISTRIBUTION IS THE C DISTRIBUTION OF THE NUMBER OF EVENTS C IN THE INTERVAL (0,ALAMBA) WHEN C THE WAITING TIME BETWEEN EVENTS C IS EXPONENTIALLY DISTRIBUTED C WITH MEAN = 1 AND STANDARD DEVIATION = 1. 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 --ALAMBA = THE SINGLE PRECISION VALUE C OF THE TAIL LENGTH PARAMETER. C ALAMBA SHOULD BE POSITIVE. 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 POISSON DISTRIBUTION C WITH TAIL LENGTH PARAMETER = ALAMBA. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--ALAMBA SHOULD BE POSITIVE. C --P SHOULD BE BETWEEN 0.0 (INCLUSIVELY) C AND 1.0 (EXCLUSIVELY). C OTHER DATAPAC SUBROUTINES NEEDED--NORPPF, POICDF. C FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, DEXP. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION AND DOUBLE PRECISION. C LANGUAGE--ANSI FORTRAN. C COMMENT--THE SINGLE PRECISION TAIL LENGTH C PARAMETER ALAMBA IS NOT RESTRICTED C TO ONLY INTEGER VALUES. C ALAMBA CAN BE SET TO ANY POSITIVE REAL C VALUE--INTEGER OR NON-INTEGER. C --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 87-121, C ESPECIALLY PAGE 102, FORMULA 36.1. C --HASTINGS AND PEACOCK, STATISTICAL C DISTRIBUTIONS--A HANDBOOK FOR C STUDENTS AND PRACTITIONERS, 1975, C PAGES 108-113. 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 146-154. C --COX AND MILLER, THE THEORY OF STOCHASTIC C PROCESSES, 1965, PAGE 7. C --GENERAL ELECTRIC COMPANY, TABLES OF THE C INDIVIDUAL AND CUMULATIVE TERMS OF POISSON C DISTRIBUTION, 1962. C --OWEN, HANDBOOK OF STATISTICAL C TABLES, 1962, PAGES 259-261. 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 DLAMBA C IPR=6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(P.LT.0.0.OR.P.GE.1.0)GOTO50 IF(ALAMBA.LE.0.0)GOTO55 GOTO90 50 WRITE(IPR,1) WRITE(IPR,46)P PPF=0.0 RETURN 55 WRITE(IPR,15) WRITE(IPR,46)ALAMBA PPF=0.0 RETURN 90 CONTINUE 1 FORMAT(1H ,115H***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO THE 1 POIPPF SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****) 15 FORMAT(1H , 91H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 POIPPF SUBROUTINE IS NON-POSITIVE *****) 46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****) C C-----START POINT----------------------------------------------------- C DLAMBA=ALAMBA 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) PPF = 0 C IF(P.EQ.0.0)GOTO110 PF0=DEXP(-DLAMBA) IF(P.LE.PF0)GOTO110 GOTO190 110 PPF=0.0 RETURN 190 CONTINUE C C DETERMINE AN INITIAL APPROXIMATION TO THE POISSON C PERCENT POINT BY USE OF THE NORMAL APPROXIMATION C TO THE POISSON. C (SEE JOHNSON AND KOTZ, DISCRETE DISTRIBUTIONS, C PAGE 102, FORMULA 36.1). C AMEAN=ALAMBA SD=SQRT(ALAMBA) CALL NORPPF(P,ZPPF) X2=AMEAN-1.0+ZPPF*SD 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 POICDF(X2,ALAMBA,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 POICDF(X2,ALAMBA,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 POICDF(X2,ALAMBA,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 CCCCC 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 POISSON PROBABILITIES FOR THE C DERIVED LOWER AND UPPER BOUNDS. C X0=IX0 X1=IX1 CALL POICDF(X0,ALAMBA,P0) CALL POICDF(X1,ALAMBA,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 POICDF(X2,ALAMBA,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)ALAMBA RETURN C 222 FORMAT(1H ,43HNO UPPER BOUND FOUND AFTER 10**7 ITERATIONS) 240 FORMAT(1H ,9HIX0 = ,I8,10X,5HP0 = ,F14.7) 241 FORMAT(1H ,9HIX1 = ,I8,10X,5HP1 = ,F14.7) 242 FORMAT(1H ,9HIX2 = ,I8,10X,5HP2 = ,F14.7) 244 FORMAT(1H ,9HP = ,F14.7) 245 FORMAT(1H ,9HALAMBA = ,F14.7) 249 FORMAT(1H ,47H***** INTERNAL ERROR IN POIPPF 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