Mercurial > octave

annotate liboctave/randpoisson.c @ 11586:12df7854fa7c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
strip trailing whitespace from source files
author John W. Eaton <jwe@octave.org>
date Thu, 20 Jan 2011 17:24:59 -0500
parents fd0a3ac60b0e
children 72c96de7a403
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7019
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
1 /*
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
2
11523
fd0a3ac60b0e update copyright notices
John W. Eaton <jwe@octave.org>
parents: 10317
diff changeset
3 Copyright (C) 2006-2011 John W. Eaton
7019
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
4
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
5 This file is part of Octave.
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
6
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
7 Octave is free software; you can redistribute it and/or modify it
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
8 under the terms of the GNU General Public License as published by the
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
9 Free Software Foundation; either version 3 of the License, or (at your
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
10 option) any later version.
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
11
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
12 Octave is distributed in the hope that it will be useful, but WITHOUT
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
15 for more details.
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
16
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
17 You should have received a copy of the GNU General Public License
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
18 along with Octave; see the file COPYING. If not, see
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
19 <http://www.gnu.org/licenses/>.
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
20
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
21 */
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
22
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
23 /* Original version written by Paul Kienzle distributed as free
4270ded9ddc6 [project @ 2007-10-13 01:42:20 by jwe]
jwe
parents: 6198
diff changeset
24 software in the in the public domain. */
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
25
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
26 /* Needs the following defines:
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
27 * NAN: value to return for Not-A-Number
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
28 * RUNI: uniform generator on (0,1)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
29 * RNOR: normal generator
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
30 * LGAMMA: log gamma function
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
31 * INFINITE: function to test whether a value is infinite
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
32 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
33
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
34 #if defined (HAVE_CONFIG_H)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
35 #include <config.h>
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
36 #endif
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
37
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
38 #include <stdio.h>
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
39
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
40 #include "f77-fcn.h"
7231
2eb392d058bb [project @ 2007-11-30 18:53:29 by jwe]
jwe
parents: 7019
diff changeset
41 #include "lo-error.h"
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
42 #include "lo-ieee.h"
7231
2eb392d058bb [project @ 2007-11-30 18:53:29 by jwe]
jwe
parents: 7019
diff changeset
43 #include "lo-math.h"
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
44 #include "randmtzig.h"
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
45 #include "randpoisson.h"
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
46
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
47 #undef NAN
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
48 #define NAN octave_NaN
6096
40212bfaed1d [project @ 2006-10-26 01:57:49 by jwe]
jwe
parents: 5775
diff changeset
49 #undef INFINITE
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
50 #define INFINITE lo_ieee_isinf
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
51 #define RUNI oct_randu()
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
52 #define RNOR oct_randn()
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
53 #define LGAMMA xlgamma
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
54
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
55 F77_RET_T
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
56 F77_FUNC (dlgams, DLGAMS) (const double *, double *, double *);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
57
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
58 static double
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
59 xlgamma (double x)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
60 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
61 double result;
7601
8a939b217863 Treat negative values to lgamma and beta correctly
David Bateman <dbateman@free.fr>
parents: 7231
diff changeset
62 #ifdef HAVE_LGAMMA
8a939b217863 Treat negative values to lgamma and beta correctly
David Bateman <dbateman@free.fr>
parents: 7231
diff changeset
63 result = lgamma (x);
8a939b217863 Treat negative values to lgamma and beta correctly
David Bateman <dbateman@free.fr>
parents: 7231
diff changeset
64 #else
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
65 double sgngam;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
66
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
67 if (lo_ieee_isnan (x))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
68 result = x;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
69 else if (x <= 0 || lo_ieee_isinf (x))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
70 result = octave_Inf;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
71 else
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
72 F77_XFCN (dlgams, DLGAMS, (&x, &result, &sgngam));
7601
8a939b217863 Treat negative values to lgamma and beta correctly
David Bateman <dbateman@free.fr>
parents: 7231
diff changeset
73 #endif
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
74 return result;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
75 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
76
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
77 /* ---- pprsc.c from Stadloeber's winrand --- */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
78
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
79 /* flogfak(k) = ln(k!) */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
80 static double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
81 flogfak (double k)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
82 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
83 #define C0 9.18938533204672742e-01
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
84 #define C1 8.33333333333333333e-02
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
85 #define C3 -2.77777777777777778e-03
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
86 #define C5 7.93650793650793651e-04
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
87 #define C7 -5.95238095238095238e-04
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
88
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
89 static double logfak[30L] = {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
90 0.00000000000000000, 0.00000000000000000, 0.69314718055994531,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
91 1.79175946922805500, 3.17805383034794562, 4.78749174278204599,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
92 6.57925121201010100, 8.52516136106541430, 10.60460290274525023,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
93 12.80182748008146961, 15.10441257307551530, 17.50230784587388584,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
94 19.98721449566188615, 22.55216385312342289, 25.19122118273868150,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
95 27.89927138384089157, 30.67186010608067280, 33.50507345013688888,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
96 36.39544520803305358, 39.33988418719949404, 42.33561646075348503,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
97 45.38013889847690803, 48.47118135183522388, 51.60667556776437357,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
98 54.78472939811231919, 58.00360522298051994, 61.26170176100200198,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
99 64.55753862700633106, 67.88974313718153498, 71.25703896716800901
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
100 };
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
101
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
102 double r, rr;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
103
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
104 if (k >= 30.0)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
105 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
106 r = 1.0 / k;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
107 rr = r * r;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
108 return ((k + 0.5)*log(k) - k + C0 + r*(C1 + rr*(C3 + rr*(C5 + rr*C7))));
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
109 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
110 else
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
111 return (logfak[(int)k]);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
112 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
113
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
114
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
115 /******************************************************************
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
116 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
117 * Poisson Distribution - Patchwork Rejection/Inversion *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
118 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
119 ******************************************************************
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
120 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
121 * For parameter my < 10 Tabulated Inversion is applied. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
122 * For my >= 10 Patchwork Rejection is employed: *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
123 * The area below the histogram function f(x) is rearranged in *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
124 * its body by certain point reflections. Within a large center *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
125 * interval variates are sampled efficiently by rejection from *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
126 * uniform hats. Rectangular immediate acceptance regions speed *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
127 * up the generation. The remaining tails are covered by *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
128 * exponential functions. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
129 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
130 ******************************************************************
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
131 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
132 * FUNCTION : - pprsc samples a random number from the Poisson *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
133 * distribution with parameter my > 0. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
134 * REFERENCE : - H. Zechner (1994): Efficient sampling from *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
135 * continuous and discrete unimodal distributions, *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
136 * Doctoral Dissertation, 156 pp., Technical *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
137 * University Graz, Austria. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
138 * SUBPROGRAM : - drand(seed) ... (0,1)-Uniform generator with *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
139 * unsigned long integer *seed. *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
140 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
141 * Implemented by H. Zechner, January 1994 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
142 * Revised by F. Niederl, July 1994 *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
143 * *
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
144 ******************************************************************/
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
145
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
146 static double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
147 f (double k, double l_nu, double c_pm)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
148 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
149 return exp(k * l_nu - flogfak(k) - c_pm);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
150 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
151
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
152 static double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
153 pprsc (double my)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
154 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
155 static double my_last = -1.0;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
156 static double m, k2, k4, k1, k5;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
157 static double dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm,
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
158 f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
159 double Dk, X, Y;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
160 double Ds, U, V, W;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
161
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
162 if (my != my_last)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
163 { /* set-up */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
164 my_last = my;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
165 /* approximate deviation of reflection points k2, k4 from my - 1/2 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
166 Ds = sqrt(my + 0.25);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
167
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
168 /* mode m, reflection points k2 and k4, and points k1 and k5, */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
169 /* which delimit the centre region of h(x) */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
170 m = floor(my);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
171 k2 = ceil(my - 0.5 - Ds);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
172 k4 = floor(my - 0.5 + Ds);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
173 k1 = k2 + k2 - m + 1L;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
174 k5 = k4 + k4 - m;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
175
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
176 /* range width of the critical left and right centre region */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
177 dl = (k2 - k1);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
178 dr = (k5 - k4);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
179
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
180 /* recurrence constants r(k)=p(k)/p(k-1) at k = k1, k2, k4+1, k5+1 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
181 r1 = my / k1;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
182 r2 = my / k2;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
183 r4 = my / (k4 + 1.0);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
184 r5 = my / (k5 + 1.0);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
185
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
186 /* reciprocal values of the scale parameters of exp. tail envelope */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
187 ll = log(r1); /* expon. tail left */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
188 lr = -log(r5); /* expon. tail right*/
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
189
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
190 /* Poisson constants, necessary for computing function values f(k) */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
191 l_my = log(my);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
192 c_pm = m * l_my - flogfak(m);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
193
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
194 /* function values f(k) = p(k)/p(m) at k = k2, k4, k1, k5 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
195 f2 = f(k2, l_my, c_pm);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
196 f4 = f(k4, l_my, c_pm);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
197 f1 = f(k1, l_my, c_pm);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
198 f5 = f(k5, l_my, c_pm);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
199
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
200 /* area of the two centre and the two exponential tail regions */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
201 /* area of the two immediate acceptance regions between k2, k4 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
202 p1 = f2 * (dl + 1.0); /* immed. left */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
203 p2 = f2 * dl + p1; /* centre left */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
204 p3 = f4 * (dr + 1.0) + p2; /* immed. right */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
205 p4 = f4 * dr + p3; /* centre right */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
206 p5 = f1 / ll + p4; /* exp. tail left */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
207 p6 = f5 / lr + p5; /* exp. tail right*/
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
208 }
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
209
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
210 for (;;)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
211 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
212 /* generate uniform number U -- U(0, p6) */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
213 /* case distinction corresponding to U */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
214 if ((U = RUNI * p6) < p2)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
215 { /* centre left */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
216
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
217 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
218 R2 = [k2, m) *[0, f2), X = k2, ... m -1 */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
219 if ((V = U - p1) < 0.0) return(k2 + floor(U/f2));
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
220 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
221 R1 = [k1, k2)*[0, f1), X = k1, ... k2-1 */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
222 if ((W = V / dl) < f1 ) return(k1 + floor(V/f1));
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
223
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
224 /* computation of candidate X < k2, and its counterpart Y > k2 */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
225 /* either squeeze-acceptance of X or acceptance-rejection of Y */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
226 Dk = floor(dl * RUNI) + 1.0;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
227 if (W <= f2 - Dk * (f2 - f2/r2))
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
228 { /* quick accept of */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
229 return(k2 - Dk); /* X = k2 - Dk */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
230 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
231 if ((V = f2 + f2 - W) < 1.0)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
232 { /* quick reject of Y*/
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
233 Y = k2 + Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
234 if (V <= f2 + Dk * (1.0 - f2)/(dl + 1.0))
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
235 { /* quick accept of */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
236 return(Y); /* Y = k2 + Dk */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
237 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
238 if (V <= f(Y, l_my, c_pm)) return(Y); /* final accept of Y*/
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
239 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
240 X = k2 - Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
241 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
242 else if (U < p4)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
243 { /* centre right */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
244 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
245 R3 = [m, k4+1)*[0, f4), X = m, ... k4 */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
246 if ((V = U - p3) < 0.0) return(k4 - floor((U - p2)/f4));
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
247 /* immediate acceptance region
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
248 R4 = [k4+1, k5+1)*[0, f5) */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
249 if ((W = V / dr) < f5 ) return(k5 - floor(V/f5));
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
250
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
251 /* computation of candidate X > k4, and its counterpart Y < k4 */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
252 /* either squeeze-acceptance of X or acceptance-rejection of Y */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
253 Dk = floor(dr * RUNI) + 1.0;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
254 if (W <= f4 - Dk * (f4 - f4*r4))
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
255 { /* quick accept of */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
256 return(k4 + Dk); /* X = k4 + Dk */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
257 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
258 if ((V = f4 + f4 - W) < 1.0)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
259 { /* quick reject of Y*/
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
260 Y = k4 - Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
261 if (V <= f4 + Dk * (1.0 - f4)/ dr)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
262 { /* quick accept of */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
263 return(Y); /* Y = k4 - Dk */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
264 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
265 if (V <= f(Y, l_my, c_pm)) return(Y); /* final accept of Y*/
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
266 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
267 X = k4 + Dk;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
268 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
269 else
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
270 {
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
271 W = RUNI;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
272 if (U < p5)
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
273 { /* expon. tail left */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
274 Dk = floor(1.0 - log(W)/ll);
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
275 if ((X = k1 - Dk) < 0L) continue; /* 0 <= X <= k1 - 1 */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
276 W *= (U - p4) * ll; /* W -- U(0, h(x)) */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
277 if (W <= f1 - Dk * (f1 - f1/r1))
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
278 return(X); /* quick accept of X*/
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
279 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
280 else
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
281 { /* expon. tail right*/
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
282 Dk = floor(1.0 - log(W)/lr);
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
283 X = k5 + Dk; /* X >= k5 + 1 */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
284 W *= (U - p5) * lr; /* W -- U(0, h(x)) */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
285 if (W <= f5 - Dk * (f5 - f5*r5))
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
286 return(X); /* quick accept of X*/
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
287 }
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
288 }
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
289
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
290 /* acceptance-rejection test of candidate X from the original area */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
291 /* test, whether W <= f(k), with W = U*h(x) and U -- U(0, 1)*/
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
292 /* log f(X) = (X - m)*log(my) - log X! + log m! */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
293 if (log(W) <= X * l_my - flogfak(X) - c_pm) return(X);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
294 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
295 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
296 /* ---- pprsc.c end ------ */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
297
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
298
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
299 /* The remainder of the file is by Paul Kienzle */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
300
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
301 /* Given uniform u, find x such that CDF(L,x)==u. Return x. */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
302 static void
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
303 poisson_cdf_lookup(double lambda, double *p, size_t n)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
304 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
305 /* Table size is predicated on the maximum value of lambda
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
306 * we want to store in the table, and the maximum value of
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
307 * returned by the uniform random number generator on [0,1).
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
308 * With lambda==10 and u_max = 1 - 1/(2^32+1), we
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
309 * have poisson_pdf(lambda,36) < 1-u_max. If instead our
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
310 * generator uses more bits of mantissa or returns a value
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
311 * in the range [0,1], then for lambda==10 we need a table
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
312 * size of 46 instead. For long doubles, the table size
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
313 * will need to be longer still. */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
314 #define TABLESIZE 46
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
315 double t[TABLESIZE];
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
316
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
317 /* Precompute the table for the u up to and including 0.458.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
318 * We will almost certainly need it. */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
319 int intlambda = (int)floor(lambda);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
320 double P;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
321 int tableidx;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
322 size_t i = n;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
323
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
324 t[0] = P = exp(-lambda);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
325 for (tableidx = 1; tableidx <= intlambda; tableidx++) {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
326 P = P*lambda/(double)tableidx;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
327 t[tableidx] = t[tableidx-1] + P;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
328 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
329
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
330 while (i-- > 0) {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
331 double u = RUNI;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
332
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
333 /* If u > 0.458 we know we can jump to floor(lambda) before
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
334 * comparing (this observation is based on Stadlober's winrand
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
335 * code). For lambda >= 1, this will be a win. Lambda < 1
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
336 * is already fast, so adding an extra comparison is not a
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
337 * problem. */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
338 int k = (u > 0.458 ? intlambda : 0);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
339
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
340 /* We aren't using a for loop here because when we find the
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
341 * right k we want to jump to the next iteration of the
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
342 * outer loop, and the continue statement will only work for
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
343 * the inner loop. */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
344 nextk:
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
345 if ( u <= t[k] ) {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
346 p[i] = (double) k;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
347 continue;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
348 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
349 if (++k < tableidx) goto nextk;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
350
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
351 /* We only need high values of the table very rarely so we
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
352 * don't automatically compute the entire table. */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
353 while (tableidx < TABLESIZE) {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
354 P = P*lambda/(double)tableidx;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
355 t[tableidx] = t[tableidx-1] + P;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
356 /* Make sure we converge to 1.0 just in case u is uniform
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
357 * on [0,1] rather than [0,1). */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
358 if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
359 tableidx++;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
360 if (u <= t[tableidx-1]) break;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
361 }
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
362
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
363 /* We are assuming that the table size is big enough here.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
364 * This should be true even if RUNI is returning values in
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
365 * the range [0,1] rather than [0,1).
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
366 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
367 p[i] = (double)(tableidx-1);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
368 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
369 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
370
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
371 /* From Press, et al., Numerical Recipes */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
372 static void
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
373 poisson_rejection (double lambda, double *p, size_t n)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
374 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
375 double sq = sqrt(2.0*lambda);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
376 double alxm = log(lambda);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
377 double g = lambda*alxm - LGAMMA(lambda+1.0);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
378 size_t i;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
379
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
380 for (i = 0; i < n; i++)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
381 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
382 double y, em, t;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
383 do {
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
384 do {
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
385 y = tan(M_PI*RUNI);
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
386 em = sq * y + lambda;
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
387 } while (em < 0.0);
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
388 em = floor(em);
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
389 t = 0.9*(1.0+y*y)*exp(em*alxm-flogfak(em)-g);
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
390 } while (RUNI > t);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
391 p[i] = em;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
392 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
393 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
394
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
395 /* The cutoff of L <= 1e8 in the following two functions before using
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
396 * the normal approximation is based on:
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
397 * > L=1e8; x=floor(linspace(0,2*L,1000));
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
398 * > max(abs(normal_pdf(x,L,L)-poisson_pdf(x,L)))
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
399 * ans = 1.1376e-28
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
400 * For L=1e7, the max is around 1e-9, which is within the step size of RUNI.
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
401 * For L>1e10 the pprsc function breaks down, as I saw from the histogram
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
402 * of a large sample, so 1e8 is both small enough and large enough. */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
403
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
404 /* Generate a set of poisson numbers with the same distribution */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
405 void
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
406 oct_fill_randp (double L, octave_idx_type n, double *p)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
407 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
408 octave_idx_type i;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
409 if (L < 0.0 || INFINITE(L))
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
410 {
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
411 for (i=0; i<n; i++)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
412 p[i] = NAN;
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
413 }
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
414 else if (L <= 10.0)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
415 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
416 poisson_cdf_lookup(L, p, n);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
417 }
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
418 else if (L <= 1e8)
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
419 {
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
420 for (i=0; i<n; i++)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
421 p[i] = pprsc(L);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
422 }
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
423 else
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
424 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
425 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
426 const double sqrtL = sqrt(L);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
427 for (i = 0; i < n; i++)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
428 {
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
429 p[i] = floor(RNOR*sqrtL + L + 0.5);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
430 if (p[i] < 0.0)
10317
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
431 p[i] = 0.0; /* will probably never happen */
42d098307c30 untabify additional source files
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
432 }
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
433 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
434 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
435
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
436 /* Generate one poisson variate */
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
437 double
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
438 oct_randp (double L)
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
439 {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
440 double ret;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
441 if (L < 0.0) ret = NAN;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
442 else if (L <= 12.0) {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
443 /* From Press, et al. Numerical recipes */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
444 double g = exp(-L);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
445 int em = -1;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
446 double t = 1.0;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
447 do {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
448 ++em;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
449 t *= RUNI;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
450 } while (t > g);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
451 ret = em;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
452 } else if (L <= 1e8) {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
453 /* numerical recipes */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
454 poisson_rejection(L, &ret, 1);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
455 } else if (INFINITE(L)) {
5775
ace8d8d26933 [project @ 2006-04-24 19:13:06 by jwe]
jwe
parents: 5742
diff changeset
456 /* FIXME R uses NaN, but the normal approx. suggests that as
5742
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
457 * limit should be inf. Which is correct? */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
458 ret = NAN;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
459 } else {
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
460 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
461 ret = floor(RNOR*sqrt(L) + L + 0.5);
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
462 if (ret < 0.0) ret = 0.0; /* will probably never happen */
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
463 }
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
464 return ret;
2cd0af543e7a [project @ 2006-04-06 08:15:49 by jwe]
jwe
parents:
diff changeset
465 }