GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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