BINRAN Subroutine

public subroutine BINRAN(N, P, Npar, Iseed, X)

NAME

binran(3f) - [M_datapac:RANDOM] generate binomial random numbers

SYNOPSIS

   SUBROUTINE BINRAN(N,P,Npar,Iseed,X)

    INTEGER,intent(in)        :: N
    REAL(kind=wp),intent(in)  :: P
    INTEGER,intent(in)        :: Npar
    INTEGER,intent(inout)     :: Iseed
    REAL(kind=wp),intent(out) :: X

DESCRIPTION

BINRAN(3f) generates a random sample of size N from the binomial
distribution with 'Bernoulli probability' parameter =
P, and integer 'number of bernoulli trials' parameter = NPAR.

The binomial distribution used herein has mean = NPAR*P and standard
deviation = sqrt(NPAR*P*(1-P)).

This distribution is defined for all discrete integer X between 0
(inclusively) and NPAR (inclusively). This distribution has the
probability function

    f(X) = c(NPAR,X) * P**X * (1-P)**(NPAR-X)

Where c(NPAR,X) is the combinatorial function equaling the number of
combinations of NPAR items taken X at a time.

The binomial distribution is the distribution of the number of
successes in NPAR Bernoulli (0,1) trials where the probability of
success in a precision trial = P.

OPTIONS

INPUT ARGUMENTS

N The desired integer number of random numbers to be generated.

P The value of the ‘Bernoulli probability’ parameter for the binomial distribution. P should be between 0.0 (exclusively) and 1.0 (exclusively).

ISEED An integer iseed value. Should be set to a non-negative value to start a new sequence of values. Will be set to -1 on return to indicate the next call should continue the current random sequence walk.

NPAR The integer value of the ‘number of Bernoulli trials’ parameter. NPAR should be a positive integer.

OUTPUT ARGUMENTS

X A vector (of dimension at least N) into which the generated random sample of size N from the binomial distribution will be placed; with ‘Bernoulli probability’ parameter = P and ‘number of Bernoulli trials’ parameter = NPAR.

EXAMPLES

Sample program:

program demo_binran
use M_datapac, only : binran
implicit none
real :: x(40), P
integer :: N, Npar, Iseed
   Iseed=0
   P=0.88
   N=size(x)
   Npar=11111
   call BINRAN(N,P,Npar,Iseed,X)
   write(*,*)X
end program demo_binran

Results:

   9746.000       9795.000       9855.000       9805.000       9787.000
   9746.000       9764.000       9774.000       9767.000       9752.000
   9770.000       9784.000       9821.000       9805.000       9784.000
   9734.000       9805.000       9813.000       9792.000       9785.000
   9784.000       9815.000       9785.000       9748.000       9718.000
   9728.000       9824.000       9782.000       9776.000       9850.000
   9770.000       9821.000       9819.000       9724.000       9783.000
   9789.000       9813.000       9798.000       9747.000       9785.000

AUTHOR

The original DATAPAC library was written by James Filliben of the
Statistical Engineering Division, National Institute of Standards
and Technology.

MAINTAINER

John Urban, 2022.05.31

LICENSE

CC0-1.0

REFERENCES

  • Johnson and Kotz, Discrete Distributions, 1969, pages 50-86.
  • Hastings and Peacock, Statistical Distributions, A Handbook for Students and Practitioners, 1975, page 41.
  • Feller, An Introduction to Probability Theory and Its Applications, Volume 1, Edition 2, 1957, pages 135-142.
  • National Bureau of Standards Applied Mathematics Series 55, 1964, page 929.
  • Kendall and Stuart, The Advanced Theory of Statistics, Volume 1, Edition 2, 1963, pages 120-125.
  • Mood and Grable, Introduction to the Theory of Statistics, Edition 2, 1963, pages 64-69.
  • Tocher, The Art Of Simulation, 1963, pages 39-40.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: N
real(kind=wp), intent(in) :: P
integer, intent(in) :: Npar
integer, intent(inout) :: Iseed
real(kind=wp), intent(out) :: X(:)

Source Code

SUBROUTINE BINRAN(N,P,Npar,Iseed,X)
INTEGER,intent(in)        :: N
REAL(kind=wp),intent(in)  :: P
INTEGER,intent(in)        :: Npar
INTEGER,intent(inout)     :: Iseed
REAL(kind=wp),intent(out) :: X(:)

REAL(kind=wp) :: g(1) , u(1)
INTEGER       :: i , ig , isum , j
!
!   NOTE THAT EVEN THOUGH THE OUTPUT FROM THIS DISCRETE RANDOM NUMBER
!   GENERATOR MUST NECESSARILY BE A SEQUENCE OF ***INTEGER*** VALUES,
!   THE OUTPUT VECTOR X IS SINGLE PRECISION IN MODE.  X HAS BEEN SPECIFIED
!   AS SINGLE PRECISION SO AS TO CONFORM WITH THE DATAPAC CONVENTION THAT
!   ALL OUTPUT VECTORS FROM ALL THIS CONVENTION IS BASED ON THE BELIEF THAT
!
!    1) A MIXTURE OF MODES (FLOATING POINT VERSUS INTEGER) IS INCONSISTENT
!       AND AN UNNECESSARY COMPLICATION IN A DATA ANALYSIS; AND
!    2) FLOATING POINT MACHINE ARITHMETIC (AS OPPOSED TO INTEGER
!       ARITHMETIC) IS THE MORE NATURAL MODE FOR DOING DATA ANALYSIS.
!
!-----START POINT-----------------------------------------------------
!
!     CHECK THE INPUT ARGUMENTS FOR ERRORS
!
   IF ( N<1 ) THEN
      WRITE (G_IO,99001)
      99001 FORMAT (' ***** FATAL ERROR--THE FIRST INPUT ARGUMENT TO BINRAN(3f) IS NON-POSITIVE *****')
      WRITE (G_IO,99005) N
      RETURN
   ELSEIF ( P<=0.0_wp .OR. P>=1.0_wp ) THEN
      WRITE (G_IO,99002)
      99002 FORMAT (' ***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO BINRAN(3f) IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****')
      WRITE (G_IO,99003) P
      99003 FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',E15.8,' *****')
      RETURN
   ELSEIF ( Npar<1 ) THEN
      WRITE (G_IO,99004)
      99004    FORMAT (' ','***** FATAL ERROR--THE THIRD INPUT ARGUMENT TO BINRAN(3f) IS NON-POSITIVE *****')
      WRITE (G_IO,99005) Npar
      RETURN
   ELSEIF ( P<0.1_wp ) THEN
      !
      !     CHECK ON THE MAGNITUDE OF P,
      !     AND BRANCH TO THE FASTER
      !     GENERATION METHOD ACCORDINGLY.
      !
      !
      !     IF P IS SMALL,
      !     GENERATE N BINOMIAL NUMBERS
      !     USING THE FACT THAT THE
      !     WAITING TIME FOR 1 SUCCESS IN
      !     BERNOULLI TRIALS HAS A
      !     GEOMETRIC DISTRIBUTION.
      !
      DO i = 1 , N
         isum = 0
         j = 1
         DO
            CALL GEORAN(1,P,Iseed,g)
            ig = g(1) + 0.5_wp
            isum = isum + ig + 1
            IF ( isum>Npar ) THEN
               X(i) = j - 1
               EXIT
            ELSE
               j = j + 1
            ENDIF
         ENDDO
      ENDDO
      GOTO 99999
   ENDIF
   !
   !     IF P IS MODERATE OR LARGE,
   !     GENERATE N BINOMIAL RANDOM NUMBERS
   !     USING THE REJECTION METHOD.
   !
   DO i = 1 , N
      isum = 0
      DO j = 1 , Npar
         CALL UNIRAN(1,Iseed,u)
         IF ( u(1)<=P ) isum = isum + 1
      ENDDO
      X(i) = isum
   ENDDO
   RETURN
   99005 FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',I0,' *****')

99999 continue
END SUBROUTINE BINRAN