1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 1998 Ross Ihaka
4 * Copyright (C) 2000-2011 The R Core Team
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, a copy is available at
18 * http://www.r-project.org/Licenses/
19 *
20 * SYNOPSIS
21 *
22 * #include <Rmath.h>
23 * double rpois(double lambda)
24 *
25 * DESCRIPTION
26 *
27 * Random variates from the Poisson distribution.
28 *
29 * REFERENCE
30 *
31 * Ahrens, J.H. and Dieter, U. (1982).
32 * Computer generation of Poisson deviates
33 * from modified normal distributions.
34 * ACM Trans. Math. Software 8, 163-179.
35 */
36
37 #include "nmath.h"
38
39 #define a0 -0.5
40 #define a1 0.3333333
41 #define a2 -0.2500068
42 #define a3 0.2000118
43 #define a4 -0.1661269
44 #define a5 0.1421878
45 #define a6 -0.1384794
46 #define a7 0.1250060
47
48 #define one_7 0.1428571428571428571
49 #define one_12 0.0833333333333333333
50 #define one_24 0.0416666666666666667
51
52 #define repeat for(;;)
53
rpois(double mu)54 double rpois(double mu)
55 {
56 /* Factorial Table (0:9)! */
57 const static double fact[10] =
58 {
59 1., 1., 2., 6., 24., 120., 720., 5040., 40320., 362880.
60 };
61
62 /* These are static --- persistent between calls for same mu : */
63 static int l, m;
64
65 static double b1, b2, c, c0, c1, c2, c3;
66 static double pp[36], p0, p, q, s, d, omega;
67 static double big_l;/* integer "w/o overflow" */
68 static double muprev = 0., muprev2 = 0.;/*, muold = 0.*/
69
70 /* Local Vars [initialize some for -Wall]: */
71 double del, difmuk= 0., E= 0., fk= 0., fx, fy, g, px, py, t, u= 0., v, x;
72 double pois = -1.;
73 int k, kflag, big_mu, new_big_mu = FALSE;
74
75 if (!R_FINITE(mu) || mu < 0)
76 ML_ERR_return_NAN;
77
78 if (mu <= 0.)
79 return 0.;
80
81 big_mu = mu >= 10.;
82 if(big_mu)
83 new_big_mu = FALSE;
84
85 if (!(big_mu && mu == muprev)) {/* maybe compute new persistent par.s */
86
87 if (big_mu) {
88 new_big_mu = TRUE;
89 /* Case A. (recalculation of s,d,l because mu has changed):
90 * The poisson probabilities pk exceed the discrete normal
91 * probabilities fk whenever k >= m(mu).
92 */
93 muprev = mu;
94 s = sqrt(mu);
95 d = 6. * mu * mu;
96 big_l = floor(mu - 1.1484);
97 /* = an upper bound to m(mu) for all mu >= 10.*/
98 }
99 else { /* Small mu ( < 10) -- not using normal approx. */
100
101 /* Case B. (start new table and calculate p0 if necessary) */
102
103 /*muprev = 0.;-* such that next time, mu != muprev ..*/
104 if (mu != muprev) {
105 muprev = mu;
106 m = imax2(1, (int) mu);
107 l = 0; /* pp[] is already ok up to pp[l] */
108 q = p0 = p = exp(-mu);
109 }
110
111 repeat {
112 /* Step U. uniform sample for inversion method */
113 u = unif_rand();
114 if (u <= p0)
115 return 0.;
116
117 /* Step T. table comparison until the end pp[l] of the
118 pp-table of cumulative poisson probabilities
119 (0.458 > ~= pp[9](= 0.45792971447) for mu=10 ) */
120 if (l != 0) {
121 for (k = (u <= 0.458) ? 1 : imin2(l, m); k <= l; k++)
122 if (u <= pp[k])
123 return (double)k;
124 if (l == 35) /* u > pp[35] */
125 continue;
126 }
127 /* Step C. creation of new poisson
128 probabilities p[l..] and their cumulatives q =: pp[k] */
129 l++;
130 for (k = l; k <= 35; k++) {
131 p *= mu / k;
132 q += p;
133 pp[k] = q;
134 if (u <= q) {
135 l = k;
136 return (double)k;
137 }
138 }
139 l = 35;
140 } /* end(repeat) */
141 }/* mu < 10 */
142
143 } /* end {initialize persistent vars} */
144
145 /* Only if mu >= 10 : ----------------------- */
146
147 /* Step N. normal sample */
148 g = mu + s * norm_rand();/* norm_rand() ~ N(0,1), standard normal */
149
150 if (g >= 0.) {
151 pois = floor(g);
152 /* Step I. immediate acceptance if pois is large enough */
153 if (pois >= big_l)
154 return pois;
155 /* Step S. squeeze acceptance */
156 fk = pois;
157 difmuk = mu - fk;
158 u = unif_rand(); /* ~ U(0,1) - sample */
159 if (d * u >= difmuk * difmuk * difmuk)
160 return pois;
161 }
162
163 /* Step P. preparations for steps Q and H.
164 (recalculations of parameters if necessary) */
165
166 if (new_big_mu || mu != muprev2) {
167 /* Careful! muprev2 is not always == muprev
168 because one might have exited in step I or S
169 */
170 muprev2 = mu;
171 omega = M_1_SQRT_2PI / s;
172 /* The quantities b1, b2, c3, c2, c1, c0 are for the Hermite
173 * approximations to the discrete normal probabilities fk. */
174
175 b1 = one_24 / mu;
176 b2 = 0.3 * b1 * b1;
177 c3 = one_7 * b1 * b2;
178 c2 = b2 - 15. * c3;
179 c1 = b1 - 6. * b2 + 45. * c3;
180 c0 = 1. - b1 + 3. * b2 - 15. * c3;
181 c = 0.1069 / mu; /* guarantees majorization by the 'hat'-function. */
182 }
183
184 if (g >= 0.) {
185 /* 'Subroutine' F is called (kflag=0 for correct return) */
186 kflag = 0;
187 goto Step_F;
188 }
189
190
191 repeat {
192 /* Step E. Exponential Sample */
193
194 E = exp_rand(); /* ~ Exp(1) (standard exponential) */
195
196 /* sample t from the laplace 'hat'
197 (if t <= -0.6744 then pk < fk for all mu >= 10.) */
198 u = 2 * unif_rand() - 1.;
199 t = 1.8 + fsign(E, u);
200 if (t > -0.6744) {
201 pois = floor(mu + s * t);
202 fk = pois;
203 difmuk = mu - fk;
204
205 /* 'subroutine' F is called (kflag=1 for correct return) */
206 kflag = 1;
207
208 Step_F: /* 'subroutine' F : calculation of px,py,fx,fy. */
209
210 if (pois < 10) { /* use factorials from table fact[] */
211 px = -mu;
212 py = pow(mu, pois) / fact[(int)pois];
213 }
214 else {
215 /* Case pois >= 10 uses polynomial approximation
216 a0-a7 for accuracy when advisable */
217 del = one_12 / fk;
218 del = del * (1. - 4.8 * del * del);
219 v = difmuk / fk;
220 if (fabs(v) <= 0.25)
221 px = fk * v * v * (((((((a7 * v + a6) * v + a5) * v + a4) *
222 v + a3) * v + a2) * v + a1) * v + a0)
223 - del;
224 else /* |v| > 1/4 */
225 px = fk * log(1. + v) - difmuk - del;
226 py = M_1_SQRT_2PI / sqrt(fk);
227 }
228 x = (0.5 - difmuk) / s;
229 x *= x;/* x^2 */
230 fx = -0.5 * x;
231 fy = omega * (((c3 * x + c2) * x + c1) * x + c0);
232 if (kflag > 0) {
233 /* Step H. Hat acceptance (E is repeated on rejection) */
234 if (c * fabs(u) <= py * exp(px + E) - fy * exp(fx + E))
235 break;
236 } else
237 /* Step Q. Quotient acceptance (rare case) */
238 if (fy - u * fy <= py * exp(px - fx))
239 break;
240 }/* t > -.67.. */
241 }
242 return pois;
243 }
244