Mercurial > octave-nkf

annotate libcruft/ranlib/ignnbn.f @ 10103:0e71ead7359d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use CALL XSTOPX instead of STOP in Fortran ranlib routines New call allows Octave's error handler to intercept otherwise fatal errors in Fortran code
author Rik <rdrider0-list@yahoo.com>
date Wed, 13 Jan 2010 22:22:46 -0800
parents df7c57a6639d
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
1 INTEGER FUNCTION ignnbn(n,p)
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
2 C**********************************************************************
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
3 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
4 C INTEGER FUNCTION IGNNBN( N, P )
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
5 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
6 C GENerate Negative BiNomial random deviate
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
7 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
8 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
9 C Function
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
10 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
11 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
12 C Generates a single random deviate from a negative binomial
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
13 C distribution.
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
14 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
15 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
16 C Arguments
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
17 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
18 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
19 C N --> Required number of events.
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
20 C INTEGER N
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
21 C JJV (N > 0)
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
22 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
23 C P --> The probability of an event during a Bernoulli trial.
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
24 C REAL P
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
25 C JJV (0.0 < P < 1.0)
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
26 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
27 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
28 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
29 C Method
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
30 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
31 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
32 C Algorithm from page 480 of
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
33 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
34 C Devroye, Luc
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
35 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
36 C Non-Uniform Random Variate Generation. Springer-Verlag,
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
37 C New York, 1986.
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
38 C
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
39 C**********************************************************************
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
40 C ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
41 C .. Scalar Arguments ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
42 REAL p
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
43 INTEGER n
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
44 C ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
45 C .. Local Scalars ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
46 REAL y,a,r
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
47 C ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
48 C .. External Functions ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
49 C JJV changed to call SGAMMA directly
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
50 C REAL gengam
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
51 REAL sgamma
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
52 INTEGER ignpoi
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
53 C EXTERNAL gengam,ignpoi
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
54 EXTERNAL sgamma,ignpoi
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
55 C ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
56 C .. Intrinsic Functions ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
57 INTRINSIC real
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
58 C ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
59 C .. Executable Statements ..
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
60 C Check Arguments
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
61 C JJV changed argumnet checker to abort if N <= 0
10103
0e71ead7359d Use CALL XSTOPX instead of STOP in Fortran ranlib routines
Rik <rdrider0-list@yahoo.com>
parents: 3188
diff changeset
62 IF (n.LE.0) CALL XSTOPX ('N <= 0 in IGNNBN')
0e71ead7359d Use CALL XSTOPX instead of STOP in Fortran ranlib routines
Rik <rdrider0-list@yahoo.com>
parents: 3188
diff changeset
63 IF (p.LE.0.0) CALL XSTOPX ('P <= 0.0 in IGNNBN')
0e71ead7359d Use CALL XSTOPX instead of STOP in Fortran ranlib routines
Rik <rdrider0-list@yahoo.com>
parents: 3188
diff changeset
64 IF (p.GE.1.0) CALL XSTOPX ('P >= 1.0 in IGNNBN')
3188
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
65
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
66 C Generate Y, a random gamma (n,(1-p)/p) variable
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
67 C JJV Note: the above parametrization is consistent with Devroye,
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
68 C JJV but gamma (p/(1-p),n) is the equivalent in our code
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
69 10 r = real(n)
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
70 a = p/ (1.0-p)
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
71 C y = gengam(a,r)
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
72 y = sgamma(r)/a
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
73
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
74 C Generate a random Poisson(y) variable
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
75 ignnbn = ignpoi(y)
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
76 RETURN
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
77
df7c57a6639d [project @ 1998-10-15 06:02:21 by jwe]
jwe
parents:
diff changeset
78 END