Chapter 18: Random Number Generation

RNGDA

Generates pseudorandom numbers from a general discrete distribution using an alias method.

Required Arguments

IOPT — Indicator of whether the alias vectors are to be initialized.   (Input)

IOPT   Action

0          The alias vectors are to be initialized using the probabilities in PROBS. IOPT is set to 0 on the first call to RNGDA.

1          The alias vectors IWK and WK are used but PROBS is not used.

IMIN — Smallest value the random deviate can assume.   (Input)
This is the value corresponding to the probability in PROBS(1).

PROBS — Vector of length NMASS containing probabilities associated with the individual mass points.   (Input)
The elements of PROBS must be nonnegative and must sum to 1.0.

IWK — Index vector of length NMASS.   (Input, if IOPT = 1; output, if IOPT = 0) IWK is a work vector.

WK — Index vector of length NMASS.   (Input, if IOPT = 1; output, if IOPT = 0) WK is a work vector.

IR — Vector of length NR containing the random discrete deviates.   (Output)

Optional Arguments

NR — Number of random numbers to generate.   (Input)
Default: NR = size (IR,1).

NMASS — Number of mass points in the discrete distribution.   (Input)
Default: NMASS = size (PROBS,1).

FORTRAN 90 Interface

Generic:                              CALL RNGDA (IOPT, IMIN, PROBS, IWK, WK, IR [,…])

Specific:                             The specific interface names are S_RNGDA and D_RNGDA.

FORTRAN 77 Interface

Single:                                CALL RNGDA (NR, IOPT, IMIN, NMASS, PROBS, IWK, WK, IR)

Double:                              The double precision name is DRNGDA.

Description

Routine RNGDA generates pseudorandom numbers from a discrete distribution with probability function given in the vector PROBS; that is

Pr(X = i) = pj

for

The algorithm is the alias method, due to Walker (1974), with modifications suggested by Kronmal and Peterson (1979). The method involves a setup phase, in which the vectors IWK and WK are filled. After the vectors are filled, the generation phase is very fast.

Comments

1.         In the interest of efficiency, this routine does only limited error checking when
IOPT = 1.

2.         The routine RNSET can be used to initialize the seed of the random number generator. The routine RNOPT can be used to select the form of the generator.

Example 1

In this example, RNGDA is used to generate five pseudorandom variates from the discrete distribution:

Pr(X = 1) = .05

Pr(X = 2) = .45

Pr(X = 3) = .31

Pr(X = 4) = .04

Pr(X = 5) = .15

When RNGDA is called the first time, IOPT is input as 0. This causes the work arrays to be initialized. In the next call, IOPT is 1, so the setup phase is bypassed.

 

      USE RNGDA_INT

      USE UMACH_INT

      USE RNSET_INT

 

      IMPLICIT   NONE

      INTEGER    NMASS, NR

      PARAMETER  (NMASS=5, NR=5)

!

      INTEGER    IMIN, IOPT, IR(NR), ISEED, IWK(NMASS), NOUT

      REAL       PROBS(NMASS), WK(NMASS)

!

      CALL UMACH (2, NOUT)

      IMIN     = 1

      PROBS(1) = 0.05

      PROBS(2) = 0.45

      PROBS(3) = 0.31

      PROBS(4) = 0.04

      PROBS(5) = 0.15

      IOPT     = 0

      ISEED    = 123457

      CALL RNSET (ISEED)

      CALL RNGDA (IOPT, IMIN, PROBS, IWK, WK, IR)

      WRITE (NOUT,99998) IR

99998 FORMAT ('          Random deviates: ', 5I4)

      IOPT = 1

      CALL RNGDA (IOPT, IMIN, PROBS, IWK, WK, IR)

      WRITE (NOUT,99999) IR

99999 FORMAT ('                           ', 5I4)

      END

Output

 

Random deviates:    3   2   2   3   5

                    1   3   4   5   3

Additional Example

Example 2

In this example, RNGDA is used to generate five pseudorandom binomial variates with parameters 20 and 0.5.

 

      USE UMACH_INT

      USE RNSET_INT

      USE RNGDA_INT

      USE BINPR_INT

 

      IMPLICIT   NONE

      INTEGER    NMASS, NR

      PARAMETER  (NMASS=21, NR=5)

!

      INTEGER    IMIN, IOPT, IR(NR), ISEED, IWK(NMASS), K, N, NOUT

      REAL       P, PROBS(NMASS), WK(NMASS)

!

      CALL UMACH (2, NOUT)

      N    = 20

      P    = 0.5

      IMIN = 0

      DO 10  K=1, NMASS

         PROBS(K) = BINPR(K-1,N,P)

   10 CONTINUE

      IOPT  = 0

      ISEED = 123457

      CALL RNSET (ISEED)

      CALL RNGDA (IOPT, IMIN, PROBS, IWK, WK, IR)

 

      WRITE (NOUT,99999) IR

99999 FORMAT ('   Binomial (20, .5) deviates: ', 5I4)

      END

Output

 

Binomial (20, .5) deviates:   12  10  16  12  11



http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260