GNU Octave  4.0.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
genmul.f
Go to the documentation of this file.
1  SUBROUTINE genmul(n,p,ncat,ix)
2 C**********************************************************************
3 C
4 C SUBROUTINE GENMUL( N, P, NCAT, IX )
5 C GENerate an observation from the MULtinomial distribution
6 C
7 C
8 C Arguments
9 C
10 C
11 C N --> Number of events that will be classified into one of
12 C the categories 1..NCAT
13 C INTEGER N
14 C
15 C P --> Vector of probabilities. P(i) is the probability that
16 C an event will be classified into category i. Thus, P(i)
17 C must be [0,1]. Only the first NCAT-1 P(i) must be defined
18 C since P(NCAT) is 1.0 minus the sum of the first
19 C NCAT-1 P(i).
20 C REAL P(NCAT-1)
21 C
22 C NCAT --> Number of categories. Length of P and IX.
23 C INTEGER NCAT
24 C
25 C IX <-- Observation from multinomial distribution. All IX(i)
26 C will be nonnegative and their sum will be N.
27 C INTEGER IX(NCAT)
28 C
29 C
30 C Method
31 C
32 C
33 C Algorithm from page 559 of
34 C
35 C Devroye, Luc
36 C
37 C Non-Uniform Random Variate Generation. Springer-Verlag,
38 C New York, 1986.
39 C
40 C**********************************************************************
41 C .. Scalar Arguments ..
42  INTEGER n,ncat
43 C ..
44 C .. Array Arguments ..
45  REAL p(*)
46  INTEGER ix(*)
47 C ..
48 C .. Local Scalars ..
49  REAL prob,ptot,sum
50  INTEGER i,icat,ntot
51 C ..
52 C .. External Functions ..
53  INTEGER ignbin
54  EXTERNAL ignbin
55 C ..
56 C .. Intrinsic Functions ..
57  INTRINSIC abs
58 C ..
59 C .. Executable Statements ..
60
61 C Check Arguments
62  IF (n.LT.0) CALL xstopx('N < 0 in GENMUL')
63  IF (ncat.LE.1) CALL xstopx('NCAT <= 1 in GENMUL')
64  ptot = 0.0
65  DO 10,i = 1,ncat - 1
66  IF (p(i).LT.0.0) CALL xstopx('Some P(i) < 0 in GENMUL')
67  IF (p(i).GT.1.0) CALL xstopx('Some P(i) > 1 in GENMUL')
68  ptot = ptot + p(i)
69  10 CONTINUE
70  IF (ptot.GT.0.99999) CALL xstopx('Sum of P(i) > 1 in GENMUL')
71
72 C Initialize variables
73  ntot = n
74  sum = 1.0
75  DO 20,i = 1,ncat
76  ix(i) = 0
77  20 CONTINUE
78
79 C Generate the observation
80  DO 30,icat = 1,ncat - 1
81  prob = p(icat)/sum
82  ix(icat) = ignbin(ntot,prob)
83  ntot = ntot - ix(icat)
84  IF (ntot.LE.0) RETURN
85  sum = sum - p(icat)
86  30 CONTINUE
87  ix(ncat) = ntot
88
89 C Finished
90  RETURN
91
92  END
subroutine genmul(n, p, ncat, ix)
Definition: genmul.f:1
T abs(T x)
Definition: pr-output.cc:3062