GNU Octave: libcruft/ranlib/genmul.f Source File

Go to the documentation of this file.
00001       SUBROUTINE genmul(n,p,ncat,ix)
00002 C**********************************************************************
00003 C
00004 C            SUBROUTINE GENMUL( N, P, NCAT, IX )
00005 C     GENerate an observation from the MULtinomial distribution
00006 C
00007 C
00008 C                              Arguments
00009 C
00010 C
00011 C     N --> Number of events that will be classified into one of
00012 C           the categories 1..NCAT
00013 C                         INTEGER N
00014 C
00015 C     P --> Vector of probabilities.  P(i) is the probability that
00016 C           an event will be classified into category i.  Thus, P(i)
00017 C           must be [0,1]. Only the first NCAT-1 P(i) must be defined
00018 C           since P(NCAT) is 1.0 minus the sum of the first
00019 C           NCAT-1 P(i).
00020 C                         REAL P(NCAT-1)
00021 C
00022 C     NCAT --> Number of categories.  Length of P and IX.
00023 C                         INTEGER NCAT
00024 C
00025 C     IX <-- Observation from multinomial distribution.  All IX(i)
00026 C            will be nonnegative and their sum will be N.
00027 C                         INTEGER IX(NCAT)
00028 C
00029 C
00030 C                              Method
00031 C
00032 C
00033 C     Algorithm from page 559 of
00034 C
00035 C     Devroye, Luc
00036 C
00037 C     Non-Uniform Random Variate Generation.  Springer-Verlag,
00038 C     New York, 1986.
00039 C
00040 C**********************************************************************
00041 C     .. Scalar Arguments ..
00042       INTEGER n,ncat
00043 C     ..
00044 C     .. Array Arguments ..
00045       REAL p(*)
00046       INTEGER ix(*)
00047 C     ..
00048 C     .. Local Scalars ..
00049       REAL prob,ptot,sum
00050       INTEGER i,icat,ntot
00051 C     ..
00052 C     .. External Functions ..
00053       INTEGER ignbin
00054       EXTERNAL ignbin
00055 C     ..
00056 C     .. Intrinsic Functions ..
00057       INTRINSIC abs
00058 C     ..
00059 C     .. Executable Statements ..
00060 
00061 C     Check Arguments
00062       IF (n.LT.0) CALL XSTOPX ('N < 0 in GENMUL')
00063       IF (ncat.LE.1) CALL XSTOPX ('NCAT <= 1 in GENMUL')
00064       ptot = 0.0
00065       DO 10,i = 1,ncat - 1
00066           IF (p(i).LT.0.0) CALL XSTOPX ('Some P(i) < 0 in GENMUL')
00067           IF (p(i).GT.1.0) CALL XSTOPX ('Some P(i) > 1 in GENMUL')
00068           ptot = ptot + p(i)
00069    10 CONTINUE
00070       IF (ptot.GT.0.99999) CALL XSTOPX ('Sum of P(i) > 1 in GENMUL')
00071 
00072 C     Initialize variables
00073       ntot = n
00074       sum = 1.0
00075       DO 20,i = 1,ncat
00076           ix(i) = 0
00077    20 CONTINUE
00078 
00079 C     Generate the observation
00080       DO 30,icat = 1,ncat - 1
00081           prob = p(icat)/sum
00082           ix(icat) = ignbin(ntot,prob)
00083           ntot = ntot - ix(icat)
00084           IF (ntot.LE.0) RETURN
00085           sum = sum - p(icat)
00086    30 CONTINUE
00087       ix(ncat) = ntot
00088 
00089 C     Finished
00090       RETURN
00091 
00092       END