| /* |
| * Mathlib : A C Library of Special Functions |
| * Copyright (C) 1998 Ross Ihaka |
| * Copyright (C) 2000-2011 The R Core Team |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or |
| * (at your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| * GNU General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, a copy is available at |
| * https://www.R-project.org/Licenses/ |
| * |
| * SYNOPSIS |
| * |
| * #include <Rmath.h> |
| * double rpois(double lambda) |
| * |
| * DESCRIPTION |
| * |
| * Random variates from the Poisson distribution. |
| * |
| * REFERENCE |
| * |
| * Ahrens, J.H. and Dieter, U. (1982). |
| * Computer generation of Poisson deviates |
| * from modified normal distributions. |
| * ACM Trans. Math. Software 8, 163-179. |
| */ |
| |
| #include "nmath.h" |
| |
| #define a0 -0.5 |
| #define a1 0.3333333 |
| #define a2 -0.2500068 |
| #define a3 0.2000118 |
| #define a4 -0.1661269 |
| #define a5 0.1421878 |
| #define a6 -0.1384794 |
| #define a7 0.1250060 |
| |
| #define one_7 0.1428571428571428571 |
| #define one_12 0.0833333333333333333 |
| #define one_24 0.0416666666666666667 |
| |
| #define repeat for(;;) |
| |
| double rpois(double mu) |
| { |
| /* Factorial Table (0:9)! */ |
| const static double fact[10] = |
| { |
| 1., 1., 2., 6., 24., 120., 720., 5040., 40320., 362880. |
| }; |
| |
| /* These are static --- persistent between calls for same mu : */ |
| static int l, m; |
| |
| static double b1, b2, c, c0, c1, c2, c3; |
| static double pp[36], p0, p, q, s, d, omega; |
| static double big_l;/* integer "w/o overflow" */ |
| static double muprev = 0., muprev2 = 0.;/*, muold = 0.*/ |
| |
| /* Local Vars [initialize some for -Wall]: */ |
| double del, difmuk= 0., E= 0., fk= 0., fx, fy, g, px, py, t, u= 0., v, x; |
| double pois = -1.; |
| int k, kflag, big_mu, new_big_mu = FALSE; |
| |
| if (!R_FINITE(mu) || mu < 0) |
| ML_ERR_return_NAN; |
| |
| if (mu <= 0.) |
| return 0.; |
| |
| big_mu = mu >= 10.; |
| if(big_mu) |
| new_big_mu = FALSE; |
| |
| if (!(big_mu && mu == muprev)) {/* maybe compute new persistent par.s */ |
| |
| if (big_mu) { |
| new_big_mu = TRUE; |
| /* Case A. (recalculation of s,d,l because mu has changed): |
| * The poisson probabilities pk exceed the discrete normal |
| * probabilities fk whenever k >= m(mu). |
| */ |
| muprev = mu; |
| s = sqrt(mu); |
| d = 6. * mu * mu; |
| big_l = floor(mu - 1.1484); |
| /* = an upper bound to m(mu) for all mu >= 10.*/ |
| } |
| else { /* Small mu ( < 10) -- not using normal approx. */ |
| |
| /* Case B. (start new table and calculate p0 if necessary) */ |
| |
| /*muprev = 0.;-* such that next time, mu != muprev ..*/ |
| if (mu != muprev) { |
| muprev = mu; |
| m = imax2(1, (int) mu); |
| l = 0; /* pp[] is already ok up to pp[l] */ |
| q = p0 = p = exp(-mu); |
| } |
| |
| repeat { |
| /* Step U. uniform sample for inversion method */ |
| u = unif_rand(); |
| if (u <= p0) |
| return 0.; |
| |
| /* Step T. table comparison until the end pp[l] of the |
| pp-table of cumulative poisson probabilities |
| (0.458 > ~= pp[9](= 0.45792971447) for mu=10 ) */ |
| if (l != 0) { |
| for (k = (u <= 0.458) ? 1 : imin2(l, m); k <= l; k++) |
| if (u <= pp[k]) |
| return (double)k; |
| if (l == 35) /* u > pp[35] */ |
| continue; |
| } |
| /* Step C. creation of new poisson |
| probabilities p[l..] and their cumulatives q =: pp[k] */ |
| l++; |
| for (k = l; k <= 35; k++) { |
| p *= mu / k; |
| q += p; |
| pp[k] = q; |
| if (u <= q) { |
| l = k; |
| return (double)k; |
| } |
| } |
| l = 35; |
| } /* end(repeat) */ |
| }/* mu < 10 */ |
| |
| } /* end {initialize persistent vars} */ |
| |
| /* Only if mu >= 10 : ----------------------- */ |
| |
| /* Step N. normal sample */ |
| g = mu + s * norm_rand();/* norm_rand() ~ N(0,1), standard normal */ |
| |
| if (g >= 0.) { |
| pois = floor(g); |
| /* Step I. immediate acceptance if pois is large enough */ |
| if (pois >= big_l) |
| return pois; |
| /* Step S. squeeze acceptance */ |
| fk = pois; |
| difmuk = mu - fk; |
| u = unif_rand(); /* ~ U(0,1) - sample */ |
| if (d * u >= difmuk * difmuk * difmuk) |
| return pois; |
| } |
| |
| /* Step P. preparations for steps Q and H. |
| (recalculations of parameters if necessary) */ |
| |
| if (new_big_mu || mu != muprev2) { |
| /* Careful! muprev2 is not always == muprev |
| because one might have exited in step I or S |
| */ |
| muprev2 = mu; |
| omega = M_1_SQRT_2PI / s; |
| /* The quantities b1, b2, c3, c2, c1, c0 are for the Hermite |
| * approximations to the discrete normal probabilities fk. */ |
| |
| b1 = one_24 / mu; |
| b2 = 0.3 * b1 * b1; |
| c3 = one_7 * b1 * b2; |
| c2 = b2 - 15. * c3; |
| c1 = b1 - 6. * b2 + 45. * c3; |
| c0 = 1. - b1 + 3. * b2 - 15. * c3; |
| c = 0.1069 / mu; /* guarantees majorization by the 'hat'-function. */ |
| } |
| |
| if (g >= 0.) { |
| /* 'Subroutine' F is called (kflag=0 for correct return) */ |
| kflag = 0; |
| goto Step_F; |
| } |
| |
| |
| repeat { |
| /* Step E. Exponential Sample */ |
| |
| E = exp_rand(); /* ~ Exp(1) (standard exponential) */ |
| |
| /* sample t from the laplace 'hat' |
| (if t <= -0.6744 then pk < fk for all mu >= 10.) */ |
| u = 2 * unif_rand() - 1.; |
| t = 1.8 + fsign(E, u); |
| if (t > -0.6744) { |
| pois = floor(mu + s * t); |
| fk = pois; |
| difmuk = mu - fk; |
| |
| /* 'subroutine' F is called (kflag=1 for correct return) */ |
| kflag = 1; |
| |
| Step_F: /* 'subroutine' F : calculation of px,py,fx,fy. */ |
| |
| if (pois < 10) { /* use factorials from table fact[] */ |
| px = -mu; |
| py = pow(mu, pois) / fact[(int)pois]; |
| } |
| else { |
| /* Case pois >= 10 uses polynomial approximation |
| a0-a7 for accuracy when advisable */ |
| del = one_12 / fk; |
| del = del * (1. - 4.8 * del * del); |
| v = difmuk / fk; |
| if (fabs(v) <= 0.25) |
| px = fk * v * v * (((((((a7 * v + a6) * v + a5) * v + a4) * |
| v + a3) * v + a2) * v + a1) * v + a0) |
| - del; |
| else /* |v| > 1/4 */ |
| px = fk * log(1. + v) - difmuk - del; |
| py = M_1_SQRT_2PI / sqrt(fk); |
| } |
| x = (0.5 - difmuk) / s; |
| x *= x;/* x^2 */ |
| fx = -0.5 * x; |
| fy = omega * (((c3 * x + c2) * x + c1) * x + c0); |
| if (kflag > 0) { |
| /* Step H. Hat acceptance (E is repeated on rejection) */ |
| if (c * fabs(u) <= py * exp(px + E) - fy * exp(fx + E)) |
| break; |
| } else |
| /* Step Q. Quotient acceptance (rare case) */ |
| if (fy - u * fy <= py * exp(px - fx)) |
| break; |
| }/* t > -.67.. */ |
| } |
| return pois; |
| } |