1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
|
#include "internal.h"
/*
* Ahrens, J. H., & Dieter, U. (1982).
* Computer Generation of Poisson Deviates from Modified Normal Distributions.
*/
static double factorial[10] = {1., 1., 2., 6., 24., 120., 720., 5040., 40320., 362880.};
static double coeffs[9] = {
-.500000000, +.333333333, -.249999856,
+.200011780, -.166684875, +.142187833,
-.124196313, +.125005956, -.114265030,
};
static inline
double
log1pmx(double x, double off)
{
int i;
double r, t;
if(-0.25 < x && x < 0.25) {
r = 0;
t = 1;
for(i=0;i<arrlen(coeffs);i++) {
r += coeffs[i]*t;
t *= x;
}
return x*x*r;
}
return log(1+x) - off;
}
static inline
double
procf(double mu, double s, int64 K, double *px, double *py, double *fx, double *fy)
{
double d, V, X;
double w, b1, b2, c1, c2, c3, c0, c;
w = 0.3989422804014327/s;
b1 = 0.041666666666666664/mu;
b2 = 0.3*b1*b1;
c3 = 0.14285714285714285*b1*b2;
c2 = b2 - 15.*c3;
c1 = b1 - 6.*b2 + 45.*c3;
c0 = 1 - b1 + 3.*b2 - 15.*c3;
c = .1069/mu;
if(K < 10) {
*px = -mu;
*py = pow(mu,K) / factorial[K];
}else{
d = 0.08333333333333333/K;
d = d - 4.8*d*d*d;
V = (mu - K) / K;
*px = K*log1pmx(V,mu-K) - d;
*py = 0.3989422804014327/sqrt(K);
}
X = (K - mu + 0.5)/s;
*fx = -0.5*X*X;
*fy = w*(((c3*X*X + c2)*X*X + c1)*X*X + c0);
return c;
}
static inline
uint64
bigpoisson(double mu)
{
int64 L,K;
double G,s,d,U,E,T;
double px,py,fx,fy,c;
s = sqrt(mu);
d = 6*mu*mu;
L = floor(mu - 1.1484);
stepN:
G = mu + s*rng·normal();
K = floor(G);
if(K<0)
goto stepP;
stepI:
if(K>=L)
return K;
stepS:
U = rng·random();
if(d*U >= (mu-K)*(mu-K)*(mu-K))
return K;
stepP:
if(G < 0)
goto stepE;
stepQ:
c = procf(mu, s, K, &px, &py, &fx, &fy);
stepE:
E = rng·exponential(1.0);
U = rng·random();
U = U + U - 1;
T = 1.8 + copysign(E,U);
if(T < 0.6744)
goto stepE;
K = floor(mu + s*T);
c = procf(mu, s, K, &px, &py, &fx, &fy);
stepH:
if(c*fabs(U) > (py*exp(px + E) - fy*exp(fx + E)))
goto stepE;
return K;
}
uint64
rng·poisson(double mean)
{
int64 n;
double z;
if(mean<10.0) {
for(n=0, z=rng·exponential(1.0); z<mean; ++n, z+=rng·exponential(1.0))
;
return n;
}
return bigpoisson(mean);
}
|
