summaryrefslogtreecommitdiff
path: root/modules/randlib/src/fortran/genmul.f
blob: 7b2be05c3c7a626f00d6717a0e5c582d281f2b46 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
      SUBROUTINE genmul(n,p,ncat,ix)
C**********************************************************************
C
C            SUBROUTINE GENMUL( N, P, NCAT, IX )
C     GENerate an observation from the MULtinomial distribution
C
C
C                              Arguments
C
C
C     N --> Number of events that will be classified into one of
C           the categories 1..NCAT
C                         INTEGER N
C
C     P --> Vector of probabilities.  P(i) is the probability that
C           an event will be classified into category i.  Thus, P(i)
C           must be [0,1]. Only the first NCAT-1 P(i) must be defined
C           since P(NCAT) is 1.0 minus the sum of the first
C           NCAT-1 P(i).
C                         DOUBLE PRECISION P(NCAT-1)
C
C     NCAT --> Number of categories.  Length of P and IX.
C                         INTEGER NCAT
C
C     IX <-- Observation from multinomial distribution.  All IX(i)
C            will be nonnegative and their sum will be N.
C                         INTEGER IX(NCAT)
C
C
C                              Method
C
C
C     Algorithm from page 559 of
C
C     Devroye, Luc
C
C     Non-Uniform Random Variate Generation.  Springer-Verlag,
C     New York, 1986.
C
C**********************************************************************
C     .. Scalar Arguments ..
      INTEGER n,ncat
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION p(*)
      INTEGER ix(*)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION prob,sum
      INTEGER i,icat,ntot
C     ..
C     .. External Functions ..
      INTEGER ignbin
      EXTERNAL ignbin
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs
C     ..
C     .. Executable Statements ..

C     Check Arguments
C     see Rand.c 
C     Initialize variables
      ntot = n
      sum = 1.0
      DO 20,i = 1,ncat
          ix(i) = 0
   20 CONTINUE

C     Generate the observation
      DO 30,icat = 1,ncat - 1
          prob = p(icat)/sum
          ix(icat) = ignbin(ntot,prob)
          ntot = ntot - ix(icat)
          IF (ntot.LE.0) RETURN
          sum = sum - p(icat)
   30 CONTINUE
      ix(ncat) = ntot

C     Finished
      RETURN

      END