1 /*
2 Copyright � 1999 CERN - European Organization for Nuclear Research.
3 Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
4 is hereby granted without fee, provided that the above copyright notice appear in all copies and
5 that both that copyright notice and this permission notice appear in supporting documentation.
6 CERN makes no representations about the suitability of this software for any purpose.
7 It is provided "as is" without expressed or implied warranty.
8 */
9 package cern.jet.random;
10
11 import cern.jet.math.Arithmetic;
12 import cern.jet.random.engine.RandomEngine;
13 import cern.jet.stat.Probability;
14 /**
15 * Poisson distribution (quick); See the <A HREF="http://www.cern.ch/RD11/rkb/AN16pp/node208.html#SECTION0002080000000000000000"> math definition</A>
16 * and <A HREF="http://www.statsoft.com/textbook/glosp.html#Poisson Distribution"> animated definition</A>.
17 * <p>
18 * <tt>p(k) = (mean^k / k!) * exp(-mean)</tt> for <tt>k >= 0</tt>.
19 * <p>
20 * Valid parameter ranges: <tt>mean > 0</tt>.
21 * Note: if <tt>mean <= 0.0</tt> then always returns zero.
22 * <p>
23 * Instance methods operate on a user supplied uniform random number generator; they are unsynchronized.
24 * <dt>
25 * Static methods operate on a default uniform random number generator; they are synchronized.
26 * <p>
27 * <b>Implementation:</b> High performance implementation.
28 * Patchwork Rejection/Inversion method.
29 * <dt>This is a port of <tt>pprsc.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
30 * C-RAND's implementation, in turn, is based upon
31 * <p>
32 * H. Zechner (1994): Efficient sampling from continuous and discrete unimodal distributions,
33 * Doctoral Dissertation, 156 pp., Technical University Graz, Austria.
34 * <p>
35 * Also see
36 * <p>
37 * Stadlober E., H. Zechner (1999), <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">The patchwork rejection method for sampling from unimodal distributions</A>,
38 * to appear in ACM Transactions on Modelling and Simulation.
39 *
40 * @author wolfgang.hoschek@cern.ch
41 * @version 1.0, 09/24/99
42 */
43 public class Poisson extends AbstractDiscreteDistribution {
44 protected double mean;
45
46 // precomputed and cached values (for performance only)
47 // cache for < SWITCH_MEAN
48 protected double my_old = -1.0;
49 protected double p,q,p0;
50 protected double[] pp = new double[36];
51 protected int llll;
52
53 // cache for >= SWITCH_MEAN
54 protected double my_last = -1.0;
55 protected double ll;
56 protected int k2, k4, k1, k5;
57 protected double dl, dr, r1, r2, r4, r5, lr, l_my, c_pm;
58 protected double f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
59
60 // cache for both;
61 protected int m;
62
63
64 protected static final double MEAN_MAX = Integer.MAX_VALUE; // for all means larger than that, we don't try to compute a poisson deviation, but return the mean.
65 protected static final double SWITCH_MEAN = 10.0; // switch from method A to method B
66
67
68 // The uniform random number generated shared by all <b>static</b> methods.
69 protected static Poisson shared = new Poisson(0.0,makeDefaultGenerator());
70 /**
71 * Constructs a poisson distribution.
72 * Example: mean=1.0.
73 */
Poisson(double mean, RandomEngine randomGenerator)74 public Poisson(double mean, RandomEngine randomGenerator) {
75 setRandomGenerator(randomGenerator);
76 setMean(mean);
77 }
78 /**
79 * Returns the cumulative distribution function.
80 */
cdf(int k)81 public double cdf(int k) {
82 return Probability.poisson(k,this.mean);
83 }
84 /**
85 * Returns a deep copy of the receiver; the copy will produce identical sequences.
86 * After this call has returned, the copy and the receiver have equal but separate state.
87 *
88 * @return a copy of the receiver.
89 */
clone()90 public Object clone() {
91 Poisson copy = (Poisson) super.clone();
92 if (this.pp != null) copy.pp = (double[]) this.pp.clone();
93 return copy;
94 }
f(int k, double l_nu, double c_pm)95 private static double f(int k, double l_nu, double c_pm) {
96 return Math.exp(k * l_nu - Arithmetic.logFactorial(k) - c_pm);
97 }
98 /**
99 * Returns a random number from the distribution.
100 */
nextInt()101 public int nextInt() {
102 return nextInt(this.mean);
103 }
104 /**
105 * Returns a random number from the distribution; bypasses the internal state.
106 */
nextInt(double theMean)107 public int nextInt(double theMean) {
108 /******************************************************************
109 * *
110 * Poisson Distribution - Patchwork Rejection/Inversion *
111 * *
112 ******************************************************************
113 * *
114 * For parameter my < 10 Tabulated Inversion is applied. *
115 * For my >= 10 Patchwork Rejection is employed: *
116 * The area below the histogram function f(x) is rearranged in *
117 * its body by certain point reflections. Within a large center *
118 * interval variates are sampled efficiently by rejection from *
119 * uniform hats. Rectangular immediate acceptance regions speed *
120 * up the generation. The remaining tails are covered by *
121 * exponential functions. *
122 * *
123 *****************************************************************/
124 RandomEngine gen = this.randomGenerator;
125 double my = theMean;
126
127 double t,g,my_k;
128
129 double gx,gy,px,py,e,x,xx,delta,v;
130 int sign;
131
132 //static double p,q,p0,pp[36];
133 //static long ll,m;
134 double u;
135 int k,i;
136
137 if (my < SWITCH_MEAN) { // CASE B: Inversion- start new table and calculate p0
138 if (my != my_old) {
139 my_old = my;
140 llll = 0;
141 p = Math.exp(-my);
142 q = p;
143 p0 = p;
144 //for (k=pp.length; --k >=0; ) pp[k] = 0;
145 }
146 m = (my > 1.0) ? (int)my : 1;
147 for(;;) {
148 u = gen.raw(); // Step U. Uniform sample
149 k = 0;
150 if (u <= p0) return(k);
151 if (llll != 0) { // Step T. Table comparison
152 i = (u > 0.458) ? Math.min(llll,m) : 1;
153 for (k = i; k <=llll; k++) if (u <= pp[k]) return(k);
154 if (llll == 35) continue;
155 }
156 for (k = llll +1; k <= 35; k++) { // Step C. Creation of new prob.
157 p *= my/(double)k;
158 q += p;
159 pp[k] = q;
160 if (u <= q) {
161 llll = k;
162 return(k);
163 }
164 }
165 llll = 35;
166 }
167 } // end my < SWITCH_MEAN
168 else if (my < MEAN_MAX ) { // CASE A: acceptance complement
169 //static double my_last = -1.0;
170 //static long int m, k2, k4, k1, k5;
171 //static double dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm,
172 // f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
173 int Dk, X, Y;
174 double Ds, U, V, W;
175
176 m = (int) my;
177 if (my != my_last) { // set-up
178 my_last = my;
179
180 // approximate deviation of reflection points k2, k4 from my - 1/2
181 Ds = Math.sqrt(my + 0.25);
182
183 // mode m, reflection points k2 and k4, and points k1 and k5, which
184 // delimit the centre region of h(x)
185 k2 = (int) Math.ceil(my - 0.5 - Ds);
186 k4 = (int) (my - 0.5 + Ds);
187 k1 = k2 + k2 - m + 1;
188 k5 = k4 + k4 - m;
189
190 // range width of the critical left and right centre region
191 dl = (double) (k2 - k1);
192 dr = (double) (k5 - k4);
193
194 // recurrence constants r(k) = p(k)/p(k-1) at k = k1, k2, k4+1, k5+1
195 r1 = my / (double) k1;
196 r2 = my / (double) k2;
197 r4 = my / (double)(k4 + 1);
198 r5 = my / (double)(k5 + 1);
199
200 // reciprocal values of the scale parameters of expon. tail envelopes
201 ll = Math.log(r1); // expon. tail left
202 lr = -Math.log(r5); // expon. tail right
203
204 // Poisson constants, necessary for computing function values f(k)
205 l_my = Math.log(my);
206 c_pm = m * l_my - Arithmetic.logFactorial(m);
207
208 // function values f(k) = p(k)/p(m) at k = k2, k4, k1, k5
209 f2 = f(k2, l_my, c_pm);
210 f4 = f(k4, l_my, c_pm);
211 f1 = f(k1, l_my, c_pm);
212 f5 = f(k5, l_my, c_pm);
213
214 // area of the two centre and the two exponential tail regions
215 // area of the two immediate acceptance regions between k2, k4
216 p1 = f2 * (dl + 1.0); // immed. left
217 p2 = f2 * dl + p1; // centre left
218 p3 = f4 * (dr + 1.0) + p2; // immed. right
219 p4 = f4 * dr + p3; // centre right
220 p5 = f1 / ll + p4; // expon. tail left
221 p6 = f5 / lr + p5; // expon. tail right
222 } // end set-up
223
224 for (;;) {
225 // generate uniform number U -- U(0, p6)
226 // case distinction corresponding to U
227 if ((U = gen.raw() * p6) < p2) { // centre left
228
229 // immediate acceptance region R2 = [k2, m) *[0, f2), X = k2, ... m -1
230 if ((V = U - p1) < 0.0) return(k2 + (int)(U/f2));
231 // immediate acceptance region R1 = [k1, k2)*[0, f1), X = k1, ... k2-1
232 if ((W = V / dl) < f1 ) return(k1 + (int)(V/f1));
233
234 // computation of candidate X < k2, and its counterpart Y > k2
235 // either squeeze-acceptance of X or acceptance-rejection of Y
236 Dk = (int)(dl * gen.raw()) + 1;
237 if (W <= f2 - Dk * (f2 - f2/r2)) { // quick accept of
238 return(k2 - Dk); // X = k2 - Dk
239 }
240 if ((V = f2 + f2 - W) < 1.0) { // quick reject of Y
241 Y = k2 + Dk;
242 if (V <= f2 + Dk * (1.0 - f2)/(dl + 1.0)) {// quick accept of
243 return(Y); // Y = k2 + Dk
244 }
245 if (V <= f(Y, l_my, c_pm)) return(Y); // final accept of Y
246 }
247 X = k2 - Dk;
248 }
249 else if (U < p4) { // centre right
250 // immediate acceptance region R3 = [m, k4+1)*[0, f4), X = m, ... k4
251 if ((V = U - p3) < 0.0) return(k4 - (int)((U - p2)/f4));
252 // immediate acceptance region R4 = [k4+1, k5+1)*[0, f5)
253 if ((W = V / dr) < f5 ) return(k5 - (int)(V/f5));
254
255 // computation of candidate X > k4, and its counterpart Y < k4
256 // either squeeze-acceptance of X or acceptance-rejection of Y
257 Dk = (int)(dr * gen.raw()) + 1;
258 if (W <= f4 - Dk * (f4 - f4*r4)) { // quick accept of
259 return(k4 + Dk); // X = k4 + Dk
260 }
261 if ((V = f4 + f4 - W) < 1.0) { // quick reject of Y
262 Y = k4 - Dk;
263 if (V <= f4 + Dk * (1.0 - f4)/ dr) { // quick accept of
264 return(Y); // Y = k4 - Dk
265 }
266 if (V <= f(Y, l_my, c_pm)) return(Y); // final accept of Y
267 }
268 X = k4 + Dk;
269 }
270 else {
271 W = gen.raw();
272 if (U < p5) { // expon. tail left
273 Dk = (int)(1.0 - Math.log(W)/ll);
274 if ((X = k1 - Dk) < 0) continue; // 0 <= X <= k1 - 1
275 W *= (U - p4) * ll; // W -- U(0, h(x))
276 if (W <= f1 - Dk * (f1 - f1/r1)) return(X); // quick accept of X
277 }
278 else { // expon. tail right
279 Dk = (int)(1.0 - Math.log(W)/lr);
280 X = k5 + Dk; // X >= k5 + 1
281 W *= (U - p5) * lr; // W -- U(0, h(x))
282 if (W <= f5 - Dk * (f5 - f5*r5)) return(X); // quick accept of X
283 }
284 }
285
286 // acceptance-rejection test of candidate X from the original area
287 // test, whether W <= f(k), with W = U*h(x) and U -- U(0, 1)
288 // log f(X) = (X - m)*log(my) - log X! + log m!
289 if (Math.log(W) <= X * l_my - Arithmetic.logFactorial(X) - c_pm) return(X);
290 }
291 }
292 else { // mean is too large
293 return (int) my;
294 }
295 }
296 /**
297 * Returns the probability distribution function.
298 */
pdf(int k)299 public double pdf(int k) {
300 return Math.exp(k*Math.log(this.mean) - Arithmetic.logFactorial(k) - this.mean);
301
302 // Overflow sensitive:
303 // return (Math.pow(mean,k) / cephes.Arithmetic.factorial(k)) * Math.exp(-this.mean);
304 }
305 /**
306 * Sets the mean.
307 */
setMean(double mean)308 public void setMean(double mean) {
309 this.mean = mean;
310 }
311 /**
312 * Returns a random number from the distribution with the given mean.
313 */
staticNextInt(double mean)314 public static int staticNextInt(double mean) {
315 synchronized (shared) {
316 shared.setMean(mean);
317 return shared.nextInt();
318 }
319 }
320 /**
321 * Returns a String representation of the receiver.
322 */
toString()323 public String toString() {
324 return this.getClass().getName()+"("+mean+")";
325 }
326 /**
327 * Sets the uniform random number generated shared by all <b>static</b> methods.
328 * @param randomGenerator the new uniform random number generator to be shared.
329 */
xstaticSetRandomGenerator(RandomEngine randomGenerator)330 private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
331 synchronized (shared) {
332 shared.setRandomGenerator(randomGenerator);
333 }
334 }
335 }
336