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