SUBROUTINE NBRAN(N,P,NPAR,ISTART,X) C C PURPOSE--THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N C FROM THE NEGATIVE BINOMIAL DISTRIBUTION C WITH SINGLE PRECISION 'BERNOULLI PROBABILITY' C PARAMETER = P, C AND INTEGER 'NUMBER OF SUCCESSES IN BERNOULLI TRIALS' C PARAMETER = NPAR. C THE NEGATIVE BINOMIAL DISTRIBUTION USED C HEREIN HAS MEAN = NPAR*(1-P)/P C AND STANDARD DEVIATION = SQRT(NPAR*(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(NPAR+X-1,NPAR) * P**NPAR * (1-P)**X. C WHERE C(NPAR+X-1,NPAR) IS THE COMBINATORIAL FUNCTION C EQUALING THE NUMBER OF COMBINATIONS OF NPAR+X-1 ITEMS C TAKEN NPAR AT A TIME. C THE NEGATIVE BINOMIAL DISTRIBUTION IS THE C DISTRIBUTION OF THE NUMBER OF FAILURES C BEFORE OBTAINING NPAR 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--N = THE DESIRED INTEGER NUMBER C OF RANDOM NUMBERS TO BE C GENERATED. 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 --NPAR = THE INTEGER VALUE C OF THE 'NUMBER OF SUCCESSES C IN BERNOULLI TRIALS' PARAMETER. C NPAR SHOULD BE A POSITIVE INTEGER. C --ISTART = AN INTEGER FLAG CODE WHICH C (IF SET TO 0) WILL START THE C GENERATOR OVER AND HENCE C PRODUCE THE SAME RANDOM SAMPLE C OVER AND OVER AGAIN C UPON SUCCESSIVE CALLS TO C THIS SUBROUTINE WITHIN A RUN; OR C (IF SET TO SOME INTEGER C VALUE NOT EQUAL TO 0, C LIKE, SAY, 1) WILL ALLOW C THE GENERATOR TO CONTINUE C FROM WHERE IT STOPPED C AND HENCE PRODUCE DIFFERENT C RANDOM SAMPLES UPON C SUCCESSIVE CALLS TO C THIS SUBROUTINE WITHIN A RUN. C OUTPUT ARGUMENTS--X = A SINGLE PRECISION VECTOR C (OF DIMENSION AT LEAST N) C INTO WHICH THE GENERATED C RANDOM SAMPLE WILL BE PLACED. C OUTPUT--A RANDOM SAMPLE OF SIZE N C FROM THE NEGATIVE BINOMIAL DISTRIBUTION C WITH 'BERNOULLI PROBABILITY' PARAMETER = P C AND 'NUMBER OF SUCCESSES IN BERNOULLI TRIALS' C PARAMETER = NPAR. C PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. C RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE C OF N FOR THIS SUBROUTINE. C --P SHOULD BE BETWEEN 0.0 (EXCLUSIVELY) C AND 1.0 (EXCLUSIVELY). C --NPAR SHOULD BE A POSITIVE INTEGER. C OTHER DATAPAC SUBROUTINES NEEDED--UNIRAN, BINRAN, GEORAN. C FORTRAN LIBRARY SUBROUTINES NEEDED--NONE. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN. C COMMENT--NOTE THAT EVEN THOUGH THE OUTPUT C FROM THIS DISCRETE RANDOM NUMBER C GENERATOR MUST NECESSARILY BE A C SEQUENCE OF ***INTEGER*** VALUES, C THE OUTPUT VECTOR 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 OUTPUT VECTORS 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--HASTINGS AND PEACOCK, STATISTICAL C DISTRIBUTIONS--A HANDBOOK FOR C STUDENTS AND PRACTITIONERS, 1975, C PAGE 95. C --JOHNSON AND KOTZ, DISCRETE C DISTRIBUTIONS, 1969, PAGES 122-142. C --FELLER, AN INTRODUCTION TO PROBABILITY C THEORY AND ITS APPLICATIONS, VOLUME 1, C EDITION 2, 1957, PAGES 155-157, 210. C --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS C SERIES 55, 1964, PAGE 929. C --KENDALL AND STUART, THE ADVANCED THEORY OF C STATISTICS, VOLUME 1, EDITION 2, 1963, PAGES 130-131. 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 DIMENSION X(1) C IPR=6 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(N.LT.1)GOTO50 IF(P.LE.0.0.OR.P.GE.1.0)GOTO55 IF(NPAR.LT.1)GOTO60 GOTO90 50 WRITE(IPR, 5) WRITE(IPR,47)N RETURN 55 WRITE(IPR,11) WRITE(IPR,46)P RETURN 60 WRITE(IPR,25) WRITE(IPR,47)NPAR RETURN 90 CONTINUE 5 FORMAT(1H , 91H***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO THE 1 BINRAN SUBROUTINE IS NON-POSITIVE *****) 11 FORMAT(1H ,115H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 BINRAN SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****) 25 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD INPUT ARGUMENT TO THE 1 BINRAN 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 CALL UNIRAN(1,ISTART,G) C C CHECK ON THE MAGNITUDE OF P, C AND BRANCH TO THE FASTER C GENERATION METHOD ACCORDINGLY. C IF(P.LT.0.1)GOTO450 C C IF P IS MODERATE OR LARGE, C GENERATE N NEGATIVE BINOMIAL NUMBERS C USING THE FACT THAT THE C WAITING TIME FOR NPAR SUCCESSES IN C BERNOULLI TRIALS HAS A C NEGATIVE BINOMIAL DISTRIBUTION. C DO100I=1,N ISUM=0 J=1 150 CALL BINRAN(1,P,1,1,B) IB=B+0.5 ISUM=ISUM+IB IF(ISUM.EQ.NPAR)GOTO250 J=J+1 GOTO150 250 X(I)=J 100 CONTINUE RETURN C C IF P IS SMALL, C GENERATE N NEGATIVE BINOMIAL NUMBERS C BY USING THE FACT THAT THE SUM C OF GEOMETRIC VARIATES IS A C NEGATIVE BINOMIAL VARIATE. C 450 DO500I=1,N ISUM=0 DO600J=1,NPAR CALL GEORAN(1,P,1,G) IG=G+0.5 ISUM=ISUM+IG 600 CONTINUE X(I)=ISUM 500 CONTINUE RETURN C END