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
ignnbn.f
Go to the documentation of this file.
1  INTEGER FUNCTION ignnbn(n,p)
2 C**********************************************************************
3 C
4 C INTEGER FUNCTION IGNNBN( N, P )
5 C
6 C GENerate Negative BiNomial random deviate
7 C
8 C
9 C Function
10 C
11 C
12 C Generates a single random deviate from a negative binomial
13 C distribution.
14 C
15 C
16 C Arguments
17 C
18 C
19 C N --> Required number of events.
20 C INTEGER N
21 C JJV (N > 0)
22 C
23 C P --> The probability of an event during a Bernoulli trial.
24 C REAL P
25 C JJV (0.0 < P < 1.0)
26 C
27 C
28 C
29 C Method
30 C
31 C
32 C Algorithm from page 480 of
33 C
34 C Devroye, Luc
35 C
36 C Non-Uniform Random Variate Generation. Springer-Verlag,
37 C New York, 1986.
38 C
39 C**********************************************************************
40 C ..
41 C .. Scalar Arguments ..
42  REAL p
43  INTEGER n
44 C ..
45 C .. Local Scalars ..
46  REAL y,a,r
47 C ..
48 C .. External Functions ..
49 C JJV changed to call SGAMMA directly
50 C REAL gengam
51  REAL sgamma
52  INTEGER ignpoi
53 C EXTERNAL gengam,ignpoi
54  EXTERNAL sgamma,ignpoi
55 C ..
56 C .. Intrinsic Functions ..
57  INTRINSIC real
58 C ..
59 C .. Executable Statements ..
60 C Check Arguments
61 C JJV changed argumnet checker to abort if N <= 0
62  IF (n.LE.0) CALL xstopx('N <= 0 in IGNNBN')
63  IF (p.LE.0.0) CALL xstopx('P <= 0.0 in IGNNBN')
64  IF (p.GE.1.0) CALL xstopx('P >= 1.0 in IGNNBN')
65 
66 C Generate Y, a random gamma (n,(1-p)/p) variable
67 C JJV Note: the above parametrization is consistent with Devroye,
68 C JJV but gamma (p/(1-p),n) is the equivalent in our code
69  10 r = real(n)
70  a = p/ (1.0-p)
71 C y = gengam(a,r)
72  y = sgamma(r)/a
73 
74 C Generate a random Poisson(y) variable
75  ignnbn = ignpoi(y)
76  RETURN
77 
78  END