File: | liboctave/numeric/randpoisson.c |
Location: | line 348, column 13 |
Description: | The right operand of '<=' is a garbage value |
1 | /* | |||
2 | ||||
3 | Copyright (C) 2006-2013 John W. Eaton | |||
4 | ||||
5 | This file is part of Octave. | |||
6 | ||||
7 | Octave is free software; you can redistribute it and/or modify it | |||
8 | under the terms of the GNU General Public License as published by the | |||
9 | Free Software Foundation; either version 3 of the License, or (at your | |||
10 | option) any later version. | |||
11 | ||||
12 | Octave is distributed in the hope that it will be useful, but WITHOUT | |||
13 | ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |||
14 | FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |||
15 | for more details. | |||
16 | ||||
17 | You should have received a copy of the GNU General Public License | |||
18 | along with Octave; see the file COPYING. If not, see | |||
19 | <http://www.gnu.org/licenses/>. | |||
20 | ||||
21 | */ | |||
22 | ||||
23 | /* Original version written by Paul Kienzle distributed as free | |||
24 | software in the in the public domain. */ | |||
25 | ||||
26 | /* Needs the following defines: | |||
27 | * NAN: value to return for Not-A-Number | |||
28 | * RUNI: uniform generator on (0,1) | |||
29 | * RNOR: normal generator | |||
30 | * LGAMMA: log gamma function | |||
31 | * INFINITE: function to test whether a value is infinite | |||
32 | */ | |||
33 | ||||
34 | #if defined (HAVE_CONFIG_H1) | |||
35 | #include <config.h> | |||
36 | #endif | |||
37 | ||||
38 | #include <stdio.h> | |||
39 | ||||
40 | #include "f77-fcn.h" | |||
41 | #include "lo-error.h" | |||
42 | #include "lo-ieee.h" | |||
43 | #include "lo-math.h" | |||
44 | #include "randmtzig.h" | |||
45 | #include "randpoisson.h" | |||
46 | ||||
47 | #undef NANoctave_NaN | |||
48 | #define NANoctave_NaN octave_NaN | |||
49 | #undef INFINITElo_ieee_isinf | |||
50 | #define INFINITElo_ieee_isinf lo_ieee_isinf | |||
51 | #define RUNIoct_randu() oct_randu() | |||
52 | #define RNORoct_randn() oct_randn() | |||
53 | #define LGAMMAxlgamma xlgamma | |||
54 | ||||
55 | F77_RET_Tint | |||
56 | F77_FUNC (dlgams, DLGAMS)dlgams_ (const double *, double *, double *); | |||
57 | ||||
58 | static double | |||
59 | xlgamma (double x) | |||
60 | { | |||
61 | double result; | |||
62 | #ifdef HAVE_LGAMMA1 | |||
63 | result = lgamma (x); | |||
64 | #else | |||
65 | double sgngam; | |||
66 | ||||
67 | if (lo_ieee_isnan (x)(sizeof (x) == sizeof (float) ? __lo_ieee_float_isnan (x) : __lo_ieee_isnan (x))) | |||
68 | result = x; | |||
69 | else if (x <= 0 || lo_ieee_isinf (x)(sizeof (x) == sizeof (float) ? __lo_ieee_float_isinf (x) : __lo_ieee_isinf (x))) | |||
70 | result = octave_Inf; | |||
71 | else | |||
72 | F77_XFCN (dlgams, DLGAMS, (&x, &result, &sgngam))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately = octave_interrupt_immediately; f77_exception_encountered = 0 ; octave_save_current_context (saved_context); if (_setjmp (current_context )) { octave_interrupt_immediately = saved_octave_interrupt_immediately ; octave_restore_current_context (saved_context); if (f77_exception_encountered ) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s" , "dlgams_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately ++; dlgams_ (&x, &result, &sgngam); octave_interrupt_immediately --; octave_restore_current_context (saved_context); } } while (0); | |||
73 | #endif | |||
74 | return result; | |||
75 | } | |||
76 | ||||
77 | /* ---- pprsc.c from Stadloeber's winrand --- */ | |||
78 | ||||
79 | /* flogfak(k) = ln(k!) */ | |||
80 | static double | |||
81 | flogfak (double k) | |||
82 | { | |||
83 | #define C09.18938533204672742e-01 9.18938533204672742e-01 | |||
84 | #define C18.33333333333333333e-02 8.33333333333333333e-02 | |||
85 | #define C3-2.77777777777777778e-03 -2.77777777777777778e-03 | |||
86 | #define C57.93650793650793651e-04 7.93650793650793651e-04 | |||
87 | #define C7-5.95238095238095238e-04 -5.95238095238095238e-04 | |||
88 | ||||
89 | static double logfak[30L] = | |||
90 | { | |||
91 | 0.00000000000000000, 0.00000000000000000, 0.69314718055994531, | |||
92 | 1.79175946922805500, 3.17805383034794562, 4.78749174278204599, | |||
93 | 6.57925121201010100, 8.52516136106541430, 10.60460290274525023, | |||
94 | 12.80182748008146961, 15.10441257307551530, 17.50230784587388584, | |||
95 | 19.98721449566188615, 22.55216385312342289, 25.19122118273868150, | |||
96 | 27.89927138384089157, 30.67186010608067280, 33.50507345013688888, | |||
97 | 36.39544520803305358, 39.33988418719949404, 42.33561646075348503, | |||
98 | 45.38013889847690803, 48.47118135183522388, 51.60667556776437357, | |||
99 | 54.78472939811231919, 58.00360522298051994, 61.26170176100200198, | |||
100 | 64.55753862700633106, 67.88974313718153498, 71.25703896716800901 | |||
101 | }; | |||
102 | ||||
103 | double r, rr; | |||
104 | ||||
105 | if (k >= 30.0) | |||
106 | { | |||
107 | r = 1.0 / k; | |||
108 | rr = r * r; | |||
109 | return ((k + 0.5)*log (k) - k + C09.18938533204672742e-01 + r*(C18.33333333333333333e-02 + rr*(C3-2.77777777777777778e-03 + rr*(C57.93650793650793651e-04 + rr*C7-5.95238095238095238e-04)))); | |||
110 | } | |||
111 | else | |||
112 | return (logfak[(int)k]); | |||
113 | } | |||
114 | ||||
115 | ||||
116 | /****************************************************************** | |||
117 | * * | |||
118 | * Poisson Distribution - Patchwork Rejection/Inversion * | |||
119 | * * | |||
120 | ****************************************************************** | |||
121 | * * | |||
122 | * For parameter my < 10 Tabulated Inversion is applied. * | |||
123 | * For my >= 10 Patchwork Rejection is employed: * | |||
124 | * The area below the histogram function f(x) is rearranged in * | |||
125 | * its body by certain point reflections. Within a large center * | |||
126 | * interval variates are sampled efficiently by rejection from * | |||
127 | * uniform hats. Rectangular immediate acceptance regions speed * | |||
128 | * up the generation. The remaining tails are covered by * | |||
129 | * exponential functions. * | |||
130 | * * | |||
131 | ****************************************************************** | |||
132 | * * | |||
133 | * FUNCTION : - pprsc samples a random number from the Poisson * | |||
134 | * distribution with parameter my > 0. * | |||
135 | * REFERENCE : - H. Zechner (1994): Efficient sampling from * | |||
136 | * continuous and discrete unimodal distributions, * | |||
137 | * Doctoral Dissertation, 156 pp., Technical * | |||
138 | * University Graz, Austria. * | |||
139 | * SUBPROGRAM : - drand(seed) ... (0,1)-Uniform generator with * | |||
140 | * unsigned long integer *seed. * | |||
141 | * * | |||
142 | * Implemented by H. Zechner, January 1994 * | |||
143 | * Revised by F. Niederl, July 1994 * | |||
144 | * * | |||
145 | ******************************************************************/ | |||
146 | ||||
147 | static double | |||
148 | f (double k, double l_nu, double c_pm) | |||
149 | { | |||
150 | return exp (k * l_nu - flogfak (k) - c_pm); | |||
151 | } | |||
152 | ||||
153 | static double | |||
154 | pprsc (double my) | |||
155 | { | |||
156 | static double my_last = -1.0; | |||
157 | static double m, k2, k4, k1, k5; | |||
158 | static double dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm, | |||
159 | f1, f2, f4, f5, p1, p2, p3, p4, p5, p6; | |||
160 | double Dk, X, Y; | |||
161 | double Ds, U, V, W; | |||
162 | ||||
163 | if (my != my_last) | |||
164 | { /* set-up */ | |||
165 | my_last = my; | |||
166 | /* approximate deviation of reflection points k2, k4 from my - 1/2 */ | |||
167 | Ds = sqrt (my + 0.25); | |||
168 | ||||
169 | /* mode m, reflection points k2 and k4, and points k1 and k5, */ | |||
170 | /* which delimit the centre region of h(x) */ | |||
171 | m = floor (my); | |||
172 | k2 = ceil (my - 0.5 - Ds); | |||
173 | k4 = floor (my - 0.5 + Ds); | |||
174 | k1 = k2 + k2 - m + 1L; | |||
175 | k5 = k4 + k4 - m; | |||
176 | ||||
177 | /* range width of the critical left and right centre region */ | |||
178 | dl = (k2 - k1); | |||
179 | dr = (k5 - k4); | |||
180 | ||||
181 | /* recurrence constants r(k)=p(k)/p(k-1) at k = k1, k2, k4+1, k5+1 */ | |||
182 | r1 = my / k1; | |||
183 | r2 = my / k2; | |||
184 | r4 = my / (k4 + 1.0); | |||
185 | r5 = my / (k5 + 1.0); | |||
186 | ||||
187 | /* reciprocal values of the scale parameters of exp. tail envelope */ | |||
188 | ll = log (r1); /* expon. tail left */ | |||
189 | lr = -log (r5); /* expon. tail right*/ | |||
190 | ||||
191 | /* Poisson constants, necessary for computing function values f(k) */ | |||
192 | l_my = log (my); | |||
193 | c_pm = m * l_my - flogfak (m); | |||
194 | ||||
195 | /* function values f(k) = p(k)/p(m) at k = k2, k4, k1, k5 */ | |||
196 | f2 = f (k2, l_my, c_pm); | |||
197 | f4 = f (k4, l_my, c_pm); | |||
198 | f1 = f (k1, l_my, c_pm); | |||
199 | f5 = f (k5, l_my, c_pm); | |||
200 | ||||
201 | /* area of the two centre and the two exponential tail regions */ | |||
202 | /* area of the two immediate acceptance regions between k2, k4 */ | |||
203 | p1 = f2 * (dl + 1.0); /* immed. left */ | |||
204 | p2 = f2 * dl + p1; /* centre left */ | |||
205 | p3 = f4 * (dr + 1.0) + p2; /* immed. right */ | |||
206 | p4 = f4 * dr + p3; /* centre right */ | |||
207 | p5 = f1 / ll + p4; /* exp. tail left */ | |||
208 | p6 = f5 / lr + p5; /* exp. tail right*/ | |||
209 | } | |||
210 | ||||
211 | for (;;) | |||
212 | { | |||
213 | /* generate uniform number U -- U(0, p6) */ | |||
214 | /* case distinction corresponding to U */ | |||
215 | if ((U = RUNIoct_randu() * p6) < p2) | |||
216 | { /* centre left */ | |||
217 | ||||
218 | /* immediate acceptance region | |||
219 | R2 = [k2, m) *[0, f2), X = k2, ... m -1 */ | |||
220 | if ((V = U - p1) < 0.0) return (k2 + floor (U/f2)); | |||
221 | /* immediate acceptance region | |||
222 | R1 = [k1, k2)*[0, f1), X = k1, ... k2-1 */ | |||
223 | if ((W = V / dl) < f1 ) return (k1 + floor (V/f1)); | |||
224 | ||||
225 | /* computation of candidate X < k2, and its counterpart Y > k2 */ | |||
226 | /* either squeeze-acceptance of X or acceptance-rejection of Y */ | |||
227 | Dk = floor (dl * RUNIoct_randu()) + 1.0; | |||
228 | if (W <= f2 - Dk * (f2 - f2/r2)) | |||
229 | { /* quick accept of */ | |||
230 | return (k2 - Dk); /* X = k2 - Dk */ | |||
231 | } | |||
232 | if ((V = f2 + f2 - W) < 1.0) | |||
233 | { /* quick reject of Y*/ | |||
234 | Y = k2 + Dk; | |||
235 | if (V <= f2 + Dk * (1.0 - f2)/(dl + 1.0)) | |||
236 | { /* quick accept of */ | |||
237 | return (Y); /* Y = k2 + Dk */ | |||
238 | } | |||
239 | if (V <= f (Y, l_my, c_pm)) return (Y); /* final accept of Y*/ | |||
240 | } | |||
241 | X = k2 - Dk; | |||
242 | } | |||
243 | else if (U < p4) | |||
244 | { /* centre right */ | |||
245 | /* immediate acceptance region | |||
246 | R3 = [m, k4+1)*[0, f4), X = m, ... k4 */ | |||
247 | if ((V = U - p3) < 0.0) return (k4 - floor ((U - p2)/f4)); | |||
248 | /* immediate acceptance region | |||
249 | R4 = [k4+1, k5+1)*[0, f5) */ | |||
250 | if ((W = V / dr) < f5 ) return (k5 - floor (V/f5)); | |||
251 | ||||
252 | /* computation of candidate X > k4, and its counterpart Y < k4 */ | |||
253 | /* either squeeze-acceptance of X or acceptance-rejection of Y */ | |||
254 | Dk = floor (dr * RUNIoct_randu()) + 1.0; | |||
255 | if (W <= f4 - Dk * (f4 - f4*r4)) | |||
256 | { /* quick accept of */ | |||
257 | return (k4 + Dk); /* X = k4 + Dk */ | |||
258 | } | |||
259 | if ((V = f4 + f4 - W) < 1.0) | |||
260 | { /* quick reject of Y*/ | |||
261 | Y = k4 - Dk; | |||
262 | if (V <= f4 + Dk * (1.0 - f4)/ dr) | |||
263 | { /* 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 | { | |||
272 | W = RUNIoct_randu(); | |||
273 | if (U < p5) | |||
274 | { /* expon. tail left */ | |||
275 | Dk = floor (1.0 - log (W)/ll); | |||
276 | if ((X = k1 - Dk) < 0L) continue; /* 0 <= X <= k1 - 1 */ | |||
277 | W *= (U - p4) * ll; /* W -- U(0, h(x)) */ | |||
278 | if (W <= f1 - Dk * (f1 - f1/r1)) | |||
279 | return (X); /* quick accept of X*/ | |||
280 | } | |||
281 | else | |||
282 | { /* expon. tail right*/ | |||
283 | Dk = floor (1.0 - log (W)/lr); | |||
284 | X = k5 + Dk; /* X >= k5 + 1 */ | |||
285 | W *= (U - p5) * lr; /* W -- U(0, h(x)) */ | |||
286 | if (W <= f5 - Dk * (f5 - f5*r5)) | |||
287 | return (X); /* quick accept of X*/ | |||
288 | } | |||
289 | } | |||
290 | ||||
291 | /* acceptance-rejection test of candidate X from the original area */ | |||
292 | /* test, whether W <= f(k), with W = U*h(x) and U -- U(0, 1)*/ | |||
293 | /* log f(X) = (X - m)*log(my) - log X! + log m! */ | |||
294 | if (log (W) <= X * l_my - flogfak (X) - c_pm) return (X); | |||
295 | } | |||
296 | } | |||
297 | /* ---- pprsc.c end ------ */ | |||
298 | ||||
299 | ||||
300 | /* The remainder of the file is by Paul Kienzle */ | |||
301 | ||||
302 | /* Given uniform u, find x such that CDF(L,x)==u. Return x. */ | |||
303 | static void | |||
304 | poisson_cdf_lookup (double lambda, double *p, size_t n) | |||
305 | { | |||
306 | /* Table size is predicated on the maximum value of lambda | |||
307 | * we want to store in the table, and the maximum value of | |||
308 | * returned by the uniform random number generator on [0,1). | |||
309 | * With lambda==10 and u_max = 1 - 1/(2^32+1), we | |||
310 | * have poisson_pdf(lambda,36) < 1-u_max. If instead our | |||
311 | * generator uses more bits of mantissa or returns a value | |||
312 | * in the range [0,1], then for lambda==10 we need a table | |||
313 | * size of 46 instead. For long doubles, the table size | |||
314 | * will need to be longer still. */ | |||
315 | #define TABLESIZE46 46 | |||
316 | double t[TABLESIZE46]; | |||
317 | ||||
318 | /* Precompute the table for the u up to and including 0.458. | |||
319 | * We will almost certainly need it. */ | |||
320 | int intlambda = (int)floor (lambda); | |||
321 | double P; | |||
322 | int tableidx; | |||
323 | size_t i = n; | |||
324 | ||||
325 | t[0] = P = exp (-lambda); | |||
326 | for (tableidx = 1; tableidx <= intlambda; tableidx++) | |||
327 | { | |||
328 | P = P*lambda/(double)tableidx; | |||
329 | t[tableidx] = t[tableidx-1] + P; | |||
330 | } | |||
331 | ||||
332 | while (i-- > 0) | |||
333 | { | |||
334 | double u = RUNIoct_randu(); | |||
335 | ||||
336 | /* If u > 0.458 we know we can jump to floor(lambda) before | |||
337 | * comparing (this observation is based on Stadlober's winrand | |||
338 | * code). For lambda >= 1, this will be a win. Lambda < 1 | |||
339 | * is already fast, so adding an extra comparison is not a | |||
340 | * problem. */ | |||
341 | int k = (u > 0.458 ? intlambda : 0); | |||
342 | ||||
343 | /* We aren't using a for loop here because when we find the | |||
344 | * right k we want to jump to the next iteration of the | |||
345 | * outer loop, and the continue statement will only work for | |||
346 | * the inner loop. */ | |||
347 | nextk: | |||
348 | if (u <= t[k]) | |||
| ||||
349 | { | |||
350 | p[i] = (double) k; | |||
351 | continue; | |||
352 | } | |||
353 | if (++k < tableidx) goto nextk; | |||
354 | ||||
355 | /* We only need high values of the table very rarely so we | |||
356 | * don't automatically compute the entire table. */ | |||
357 | while (tableidx < TABLESIZE46) | |||
358 | { | |||
359 | P = P*lambda/(double)tableidx; | |||
360 | t[tableidx] = t[tableidx-1] + P; | |||
361 | /* Make sure we converge to 1.0 just in case u is uniform | |||
362 | * on [0,1] rather than [0,1). */ | |||
363 | if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0; | |||
364 | tableidx++; | |||
365 | if (u <= t[tableidx-1]) break; | |||
366 | } | |||
367 | ||||
368 | /* We are assuming that the table size is big enough here. | |||
369 | * This should be true even if RUNI is returning values in | |||
370 | * the range [0,1] rather than [0,1). */ | |||
371 | p[i] = (double)(tableidx-1); | |||
372 | } | |||
373 | } | |||
374 | ||||
375 | static void | |||
376 | poisson_cdf_lookup_float (double lambda, float *p, size_t n) | |||
377 | { | |||
378 | double t[TABLESIZE46]; | |||
379 | ||||
380 | /* Precompute the table for the u up to and including 0.458. | |||
381 | * We will almost certainly need it. */ | |||
382 | int intlambda = (int)floor (lambda); | |||
383 | double P; | |||
384 | int tableidx; | |||
385 | size_t i = n; | |||
386 | ||||
387 | t[0] = P = exp (-lambda); | |||
388 | for (tableidx = 1; tableidx <= intlambda; tableidx++) | |||
389 | { | |||
390 | P = P*lambda/(double)tableidx; | |||
391 | t[tableidx] = t[tableidx-1] + P; | |||
392 | } | |||
393 | ||||
394 | while (i-- > 0) | |||
395 | { | |||
396 | double u = RUNIoct_randu(); | |||
397 | int k = (u > 0.458 ? intlambda : 0); | |||
398 | nextk: | |||
399 | if (u <= t[k]) | |||
400 | { | |||
401 | p[i] = (float) k; | |||
402 | continue; | |||
403 | } | |||
404 | if (++k < tableidx) goto nextk; | |||
405 | ||||
406 | while (tableidx < TABLESIZE46) | |||
407 | { | |||
408 | P = P*lambda/(double)tableidx; | |||
409 | t[tableidx] = t[tableidx-1] + P; | |||
410 | if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0; | |||
411 | tableidx++; | |||
412 | if (u <= t[tableidx-1]) break; | |||
413 | } | |||
414 | ||||
415 | p[i] = (float)(tableidx-1); | |||
416 | } | |||
417 | } | |||
418 | ||||
419 | /* From Press, et al., Numerical Recipes */ | |||
420 | static void | |||
421 | poisson_rejection (double lambda, double *p, size_t n) | |||
422 | { | |||
423 | double sq = sqrt (2.0*lambda); | |||
424 | double alxm = log (lambda); | |||
425 | double g = lambda*alxm - LGAMMAxlgamma(lambda+1.0); | |||
426 | size_t i; | |||
427 | ||||
428 | for (i = 0; i < n; i++) | |||
429 | { | |||
430 | double y, em, t; | |||
431 | do | |||
432 | { | |||
433 | do | |||
434 | { | |||
435 | y = tan (M_PI3.14159265358979323846*RUNIoct_randu()); | |||
436 | em = sq * y + lambda; | |||
437 | } while (em < 0.0); | |||
438 | em = floor (em); | |||
439 | t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g); | |||
440 | } while (RUNIoct_randu() > t); | |||
441 | p[i] = em; | |||
442 | } | |||
443 | } | |||
444 | ||||
445 | /* From Press, et al., Numerical Recipes */ | |||
446 | static void | |||
447 | poisson_rejection_float (double lambda, float *p, size_t n) | |||
448 | { | |||
449 | double sq = sqrt (2.0*lambda); | |||
450 | double alxm = log (lambda); | |||
451 | double g = lambda*alxm - LGAMMAxlgamma(lambda+1.0); | |||
452 | size_t i; | |||
453 | ||||
454 | for (i = 0; i < n; i++) | |||
455 | { | |||
456 | double y, em, t; | |||
457 | do | |||
458 | { | |||
459 | do | |||
460 | { | |||
461 | y = tan (M_PI3.14159265358979323846*RUNIoct_randu()); | |||
462 | em = sq * y + lambda; | |||
463 | } while (em < 0.0); | |||
464 | em = floor (em); | |||
465 | t = 0.9*(1.0+y*y)*exp (em*alxm-flogfak (em)-g); | |||
466 | } while (RUNIoct_randu() > t); | |||
467 | p[i] = em; | |||
468 | } | |||
469 | } | |||
470 | ||||
471 | /* The cutoff of L <= 1e8 in the following two functions before using | |||
472 | * the normal approximation is based on: | |||
473 | * > L=1e8; x=floor(linspace(0,2*L,1000)); | |||
474 | * > max(abs(normal_pdf(x,L,L)-poisson_pdf(x,L))) | |||
475 | * ans = 1.1376e-28 | |||
476 | * For L=1e7, the max is around 1e-9, which is within the step size of RUNI. | |||
477 | * For L>1e10 the pprsc function breaks down, as I saw from the histogram | |||
478 | * of a large sample, so 1e8 is both small enough and large enough. */ | |||
479 | ||||
480 | /* Generate a set of poisson numbers with the same distribution */ | |||
481 | void | |||
482 | oct_fill_randp (double L, octave_idx_type n, double *p) | |||
483 | { | |||
484 | octave_idx_type i; | |||
485 | if (L < 0.0 || INFINITE(L)(sizeof (L) == sizeof (float) ? __lo_ieee_float_isinf (L) : __lo_ieee_isinf (L))) | |||
| ||||
486 | { | |||
487 | for (i=0; i<n; i++) | |||
488 | p[i] = NANoctave_NaN; | |||
489 | } | |||
490 | else if (L <= 10.0) | |||
491 | { | |||
492 | poisson_cdf_lookup (L, p, n); | |||
493 | } | |||
494 | else if (L <= 1e8) | |||
495 | { | |||
496 | for (i=0; i<n; i++) | |||
497 | p[i] = pprsc (L); | |||
498 | } | |||
499 | else | |||
500 | { | |||
501 | /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ | |||
502 | const double sqrtL = sqrt (L); | |||
503 | for (i = 0; i < n; i++) | |||
504 | { | |||
505 | p[i] = floor (RNORoct_randn()*sqrtL + L + 0.5); | |||
506 | if (p[i] < 0.0) | |||
507 | p[i] = 0.0; /* will probably never happen */ | |||
508 | } | |||
509 | } | |||
510 | } | |||
511 | ||||
512 | /* Generate one poisson variate */ | |||
513 | double | |||
514 | oct_randp (double L) | |||
515 | { | |||
516 | double ret; | |||
517 | if (L < 0.0) ret = NANoctave_NaN; | |||
518 | else if (L <= 12.0) | |||
519 | { | |||
520 | /* From Press, et al. Numerical recipes */ | |||
521 | double g = exp (-L); | |||
522 | int em = -1; | |||
523 | double t = 1.0; | |||
524 | do | |||
525 | { | |||
526 | ++em; | |||
527 | t *= RUNIoct_randu(); | |||
528 | } while (t > g); | |||
529 | ret = em; | |||
530 | } | |||
531 | else if (L <= 1e8) | |||
532 | { | |||
533 | /* numerical recipes */ | |||
534 | poisson_rejection (L, &ret, 1); | |||
535 | } | |||
536 | else if (INFINITE(L)(sizeof (L) == sizeof (float) ? __lo_ieee_float_isinf (L) : __lo_ieee_isinf (L))) | |||
537 | { | |||
538 | /* FIXME R uses NaN, but the normal approx. suggests that as | |||
539 | * limit should be inf. Which is correct? */ | |||
540 | ret = NANoctave_NaN; | |||
541 | } | |||
542 | else | |||
543 | { | |||
544 | /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ | |||
545 | ret = floor (RNORoct_randn()*sqrt (L) + L + 0.5); | |||
546 | if (ret < 0.0) ret = 0.0; /* will probably never happen */ | |||
547 | } | |||
548 | return ret; | |||
549 | } | |||
550 | ||||
551 | /* Generate a set of poisson numbers with the same distribution */ | |||
552 | void | |||
553 | oct_fill_float_randp (float FL, octave_idx_type n, float *p) | |||
554 | { | |||
555 | double L = FL; | |||
556 | octave_idx_type i; | |||
557 | if (L < 0.0 || INFINITE(L)(sizeof (L) == sizeof (float) ? __lo_ieee_float_isinf (L) : __lo_ieee_isinf (L))) | |||
558 | { | |||
559 | for (i=0; i<n; i++) | |||
560 | p[i] = NANoctave_NaN; | |||
561 | } | |||
562 | else if (L <= 10.0) | |||
563 | { | |||
564 | poisson_cdf_lookup_float (L, p, n); | |||
565 | } | |||
566 | else if (L <= 1e8) | |||
567 | { | |||
568 | for (i=0; i<n; i++) | |||
569 | p[i] = pprsc (L); | |||
570 | } | |||
571 | else | |||
572 | { | |||
573 | /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ | |||
574 | const double sqrtL = sqrt (L); | |||
575 | for (i = 0; i < n; i++) | |||
576 | { | |||
577 | p[i] = floor (RNORoct_randn()*sqrtL + L + 0.5); | |||
578 | if (p[i] < 0.0) | |||
579 | p[i] = 0.0; /* will probably never happen */ | |||
580 | } | |||
581 | } | |||
582 | } | |||
583 | ||||
584 | /* Generate one poisson variate */ | |||
585 | float | |||
586 | oct_float_randp (float FL) | |||
587 | { | |||
588 | double L = FL; | |||
589 | float ret; | |||
590 | if (L < 0.0) ret = NANoctave_NaN; | |||
591 | else if (L <= 12.0) | |||
592 | { | |||
593 | /* From Press, et al. Numerical recipes */ | |||
594 | double g = exp (-L); | |||
595 | int em = -1; | |||
596 | double t = 1.0; | |||
597 | do | |||
598 | { | |||
599 | ++em; | |||
600 | t *= RUNIoct_randu(); | |||
601 | } while (t > g); | |||
602 | ret = em; | |||
603 | } | |||
604 | else if (L <= 1e8) | |||
605 | { | |||
606 | /* numerical recipes */ | |||
607 | poisson_rejection_float (L, &ret, 1); | |||
608 | } | |||
609 | else if (INFINITE(L)(sizeof (L) == sizeof (float) ? __lo_ieee_float_isinf (L) : __lo_ieee_isinf (L))) | |||
610 | { | |||
611 | /* FIXME R uses NaN, but the normal approx. suggests that as | |||
612 | * limit should be inf. Which is correct? */ | |||
613 | ret = NANoctave_NaN; | |||
614 | } | |||
615 | else | |||
616 | { | |||
617 | /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ | |||
618 | ret = floor (RNORoct_randn()*sqrt (L) + L + 0.5); | |||
619 | if (ret < 0.0) ret = 0.0; /* will probably never happen */ | |||
620 | } | |||
621 | return ret; | |||
622 | } |