SUBROUTINE BINRAN(N,P,NPAR,ISEED,X) C C PURPOSE--THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N C FROM THE BINOMIAL DISTRIBUTION C WITH SINGLE PRECISION 'BERNOULLI PROBABILITY' C PARAMETER = P, C AND INTEGER 'NUMBER OF BERNOULLI TRIALS' C PARAMETER = NPAR. C THE BINOMIAL DISTRIBUTION USED C HEREIN HAS MEAN = NPAR*P C AND STANDARD DEVIATION = SQRT(NPAR*P*(1-P)). C THIS DISTRIBUTION IS DEFINED FOR ALL C DISCRETE INTEGER X BETWEEN 0 (INCLUSIVELY) C AND NPAR (INCLUSIVELY). C THIS DISTRIBUTION HAS THE PROBABILITY FUNCTION C F(X) = C(NPAR,X) * P**X * (1-P)**(NPAR-X). C WHERE C(NPAR,X) IS THE COMBINATORIAL FUNCTION C EQUALING THE NUMBER OF COMBINATIONS OF NPAR ITEMS C TAKEN X AT A TIME. C THE BINOMIAL DISTRIBUTION IS THE C DISTRIBUTION OF THE NUMBER OF C SUCCESSES IN NPAR 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 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 BERNOULLI TRIALS' C PARAMETER. C NPAR SHOULD BE A POSITIVE INTEGER. 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 BINOMIAL DISTRIBUTION C WITH 'BERNOULLI PROBABILITY' PARAMETER = P C AND 'NUMBER OF BERNOULLI TRIALS' 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, GEORAN. C FORTRAN LIBRARY SUBROUTINES NEEDED--NONE. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN (1977) 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--JOHNSON AND KOTZ, DISCRETE C DISTRIBUTIONS, 1969, PAGES 50-86. C --HASTINGS AND PEACOCK, STATISTICAL C DISTRIBUTIONS--A HANDBOOK FOR C STUDENTS AND PRACTITIONERS, 1975, C PAGE 41. C --FELLER, AN INTRODUCTION TO PROBABILITY C THEORY AND ITS APPLICATIONS, VOLUME 1, C EDITION 2, 1957, PAGES 135-142. 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 120-125. C --MOOD AND GRABLE, INTRODUCTION TO THE THEORY C OF STATISTICS, EDITION 2, 1963, PAGES 64-69. C --TOCHER, THE ART OF SIMULATION, C 1963, PAGES 39-40. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING DIVISION C CENTER FOR APPLIED MATHEMATICS C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE--301-921-3651 C NOTE--DATAPLOT IS A REGISTERED TRADEMARK C OF THE NATIONAL BUREAU OF STANDARDS. C THIS SUBROUTINE MAY NOT BE COPIED, EXTRACTED, C MODIFIED, OR OTHERWISE USED IN A CONTEXT C OUTSIDE OF THE DATAPLOT LANGUAGE/SYSTEM. C LANGUAGE--ANSI FORTRAN (1966) C EXCEPTION--HOLLERITH STRINGS IN FORMAT STATEMENTS C DENOTED BY QUOTES RATHER THAN NH. C VERSION NUMBER--82/7 C ORIGINAL VERSION--NOVEMBER 1975. C UPDATED --DECEMBER 1981. C UPDATED --MAY 1982. C C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES------------------- C C--------------------------------------------------------------------- C DIMENSION X(*) C C--------------------------------------------------------------------- C CCCCC CHARACTER*4 IFEEDB CCCCC CHARACTER*4 IPRINT C CCCCC COMMON /MACH/IRD,IPR,CPUMIN,CPUMAX,NUMBPC,NUMCPW,NUMBPW CCCCC COMMON /PRINT/IFEEDB,IPRINT C IPR=6 C C-----START POINT----------------------------------------------------- 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 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 BINOMIAL RANDOM NUMBERS C USING THE REJECTION METHOD. C DO100I=1,N ISUM=0 DO200J=1,NPAR CALL UNIRAN(1,ISEED,U) IF(U.LE.P)ISUM=ISUM+1 200 CONTINUE X(I)=ISUM 100 CONTINUE RETURN C C IF P IS SMALL, C GENERATE N BINOMIAL NUMBERS C USING THE FACT THAT THE C WAITING TIME FOR 1 SUCCESS IN C BERNOULLI TRIALS HAS A C GEOMETRIC DISTRIBUTION. C 450 DO500I=1,N ISUM=0 J=1 550 CALL GEORAN(1,P,ISEED,G) IG=G+0.5 ISUM=ISUM+IG+1 IF(ISUM.GT.NPAR)GOTO650 J=J+1 GOTO550 650 X(I)=J-1 500 CONTINUE RETURN C END