Mercurial > hg > octave-nkf

comparison liboctave/randpoisson.c @ 10317:42d098307c30

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
untabify additional source files
author John W. Eaton <jwe@octave.org>
date Thu, 11 Feb 2010 13:30:42 -0500
parents 4c0cdbe0acca
children fd0a3ac60b0e
comparison
equal deleted inserted replaced
10316:9966f1f71c32 10317:42d098307c30
210 for (;;) 210 for (;;)
211 { 211 {
212 /* generate uniform number U -- U(0, p6) */ 212 /* generate uniform number U -- U(0, p6) */
213 /* case distinction corresponding to U */ 213 /* case distinction corresponding to U */
214 if ((U = RUNI * p6) < p2) 214 if ((U = RUNI * p6) < p2)
215 { /* centre left */ 215 { /* centre left */
216 216
217 /* immediate acceptance region 217 /* immediate acceptance region
218 R2 = [k2, m) *[0, f2), X = k2, ... m -1 */ 218 R2 = [k2, m) *[0, f2), X = k2, ... m -1 */
219 if ((V = U - p1) < 0.0) return(k2 + floor(U/f2)); 219 if ((V = U - p1) < 0.0) return(k2 + floor(U/f2));
220 /* immediate acceptance region 220 /* immediate acceptance region
221 R1 = [k1, k2)*[0, f1), X = k1, ... k2-1 */ 221 R1 = [k1, k2)*[0, f1), X = k1, ... k2-1 */
222 if ((W = V / dl) < f1 ) return(k1 + floor(V/f1)); 222 if ((W = V / dl) < f1 ) return(k1 + floor(V/f1));
223 223
224 /* computation of candidate X < k2, and its counterpart Y > k2 */ 224 /* computation of candidate X < k2, and its counterpart Y > k2 */
225 /* either squeeze-acceptance of X or acceptance-rejection of Y */ 225 /* either squeeze-acceptance of X or acceptance-rejection of Y */
226 Dk = floor(dl * RUNI) + 1.0; 226 Dk = floor(dl * RUNI) + 1.0;
227 if (W <= f2 - Dk * (f2 - f2/r2)) 227 if (W <= f2 - Dk * (f2 - f2/r2))
228 { /* quick accept of */ 228 { /* quick accept of */
229 return(k2 - Dk); /* X = k2 - Dk */ 229 return(k2 - Dk); /* X = k2 - Dk */
230 } 230 }
231 if ((V = f2 + f2 - W) < 1.0) 231 if ((V = f2 + f2 - W) < 1.0)
232 { /* quick reject of Y*/ 232 { /* quick reject of Y*/
233 Y = k2 + Dk; 233 Y = k2 + Dk;
234 if (V <= f2 + Dk * (1.0 - f2)/(dl + 1.0)) 234 if (V <= f2 + Dk * (1.0 - f2)/(dl + 1.0))
235 { /* quick accept of */ 235 { /* quick accept of */
236 return(Y); /* Y = k2 + Dk */ 236 return(Y); /* Y = k2 + Dk */
237 } 237 }
238 if (V <= f(Y, l_my, c_pm)) return(Y); /* final accept of Y*/ 238 if (V <= f(Y, l_my, c_pm)) return(Y); /* final accept of Y*/
239 } 239 }
240 X = k2 - Dk; 240 X = k2 - Dk;
241 } 241 }
242 else if (U < p4) 242 else if (U < p4)
243 { /* centre right */ 243 { /* centre right */
244 /* immediate acceptance region 244 /* immediate acceptance region
245 R3 = [m, k4+1)*[0, f4), X = m, ... k4 */ 245 R3 = [m, k4+1)*[0, f4), X = m, ... k4 */
246 if ((V = U - p3) < 0.0) return(k4 - floor((U - p2)/f4)); 246 if ((V = U - p3) < 0.0) return(k4 - floor((U - p2)/f4));
247 /* immediate acceptance region 247 /* immediate acceptance region
248 R4 = [k4+1, k5+1)*[0, f5) */ 248 R4 = [k4+1, k5+1)*[0, f5) */
249 if ((W = V / dr) < f5 ) return(k5 - floor(V/f5)); 249 if ((W = V / dr) < f5 ) return(k5 - floor(V/f5));
250 250
251 /* computation of candidate X > k4, and its counterpart Y < k4 */ 251 /* computation of candidate X > k4, and its counterpart Y < k4 */
252 /* either squeeze-acceptance of X or acceptance-rejection of Y */ 252 /* either squeeze-acceptance of X or acceptance-rejection of Y */
253 Dk = floor(dr * RUNI) + 1.0; 253 Dk = floor(dr * RUNI) + 1.0;
254 if (W <= f4 - Dk * (f4 - f4*r4)) 254 if (W <= f4 - Dk * (f4 - f4*r4))
255 { /* quick accept of */ 255 { /* quick accept of */
256 return(k4 + Dk); /* X = k4 + Dk */ 256 return(k4 + Dk); /* X = k4 + Dk */
257 } 257 }
258 if ((V = f4 + f4 - W) < 1.0) 258 if ((V = f4 + f4 - W) < 1.0)
259 { /* quick reject of Y*/ 259 { /* quick reject of Y*/
260 Y = k4 - Dk; 260 Y = k4 - Dk;
261 if (V <= f4 + Dk * (1.0 - f4)/ dr) 261 if (V <= f4 + Dk * (1.0 - f4)/ dr)
262 { /* quick accept of */ 262 { /* quick accept of */
263 return(Y); /* Y = k4 - Dk */ 263 return(Y); /* Y = k4 - Dk */
264 } 264 }
265 if (V <= f(Y, l_my, c_pm)) return(Y); /* final accept of Y*/ 265 if (V <= f(Y, l_my, c_pm)) return(Y); /* final accept of Y*/
266 } 266 }
267 X = k4 + Dk; 267 X = k4 + Dk;
268 } 268 }
269 else 269 else
270 { 270 {
271 W = RUNI; 271 W = RUNI;
272 if (U < p5) 272 if (U < p5)
273 { /* expon. tail left */ 273 { /* expon. tail left */
274 Dk = floor(1.0 - log(W)/ll); 274 Dk = floor(1.0 - log(W)/ll);
275 if ((X = k1 - Dk) < 0L) continue; /* 0 <= X <= k1 - 1 */ 275 if ((X = k1 - Dk) < 0L) continue; /* 0 <= X <= k1 - 1 */
276 W *= (U - p4) * ll; /* W -- U(0, h(x)) */ 276 W *= (U - p4) * ll; /* W -- U(0, h(x)) */
277 if (W <= f1 - Dk * (f1 - f1/r1)) 277 if (W <= f1 - Dk * (f1 - f1/r1))
278 return(X); /* quick accept of X*/ 278 return(X); /* quick accept of X*/
279 } 279 }
280 else 280 else
281 { /* expon. tail right*/ 281 { /* expon. tail right*/
282 Dk = floor(1.0 - log(W)/lr); 282 Dk = floor(1.0 - log(W)/lr);
283 X = k5 + Dk; /* X >= k5 + 1 */ 283 X = k5 + Dk; /* X >= k5 + 1 */
284 W *= (U - p5) * lr; /* W -- U(0, h(x)) */ 284 W *= (U - p5) * lr; /* W -- U(0, h(x)) */
285 if (W <= f5 - Dk * (f5 - f5*r5)) 285 if (W <= f5 - Dk * (f5 - f5*r5))
286 return(X); /* quick accept of X*/ 286 return(X); /* quick accept of X*/
287 } 287 }
288 } 288 }
289 289
290 /* acceptance-rejection test of candidate X from the original area */ 290 /* acceptance-rejection test of candidate X from the original area */
291 /* test, whether W <= f(k), with W = U*h(x) and U -- U(0, 1)*/ 291 /* test, whether W <= f(k), with W = U*h(x) and U -- U(0, 1)*/
292 /* log f(X) = (X - m)*log(my) - log X! + log m! */ 292 /* log f(X) = (X - m)*log(my) - log X! + log m! */
293 if (log(W) <= X * l_my - flogfak(X) - c_pm) return(X); 293 if (log(W) <= X * l_my - flogfak(X) - c_pm) return(X);
379 379
380 for (i = 0; i < n; i++) 380 for (i = 0; i < n; i++)
381 { 381 {
382 double y, em, t; 382 double y, em, t;
383 do { 383 do {
384 do { 384 do {
385 y = tan(M_PI*RUNI); 385 y = tan(M_PI*RUNI);
386 em = sq * y + lambda; 386 em = sq * y + lambda;
387 } while (em < 0.0); 387 } while (em < 0.0);
388 em = floor(em); 388 em = floor(em);
389 t = 0.9*(1.0+y*y)*exp(em*alxm-flogfak(em)-g); 389 t = 0.9*(1.0+y*y)*exp(em*alxm-flogfak(em)-g);
390 } while (RUNI > t); 390 } while (RUNI > t);
391 p[i] = em; 391 p[i] = em;
392 } 392 }
393 } 393 }
394 394
407 { 407 {
408 octave_idx_type i; 408 octave_idx_type i;
409 if (L < 0.0 || INFINITE(L)) 409 if (L < 0.0 || INFINITE(L))
410 { 410 {
411 for (i=0; i<n; i++) 411 for (i=0; i<n; i++)
412 p[i] = NAN; 412 p[i] = NAN;
413 } 413 }
414 else if (L <= 10.0) 414 else if (L <= 10.0)
415 { 415 {
416 poisson_cdf_lookup(L, p, n); 416 poisson_cdf_lookup(L, p, n);
417 } 417 }
418 else if (L <= 1e8) 418 else if (L <= 1e8)
419 { 419 {
420 for (i=0; i<n; i++) 420 for (i=0; i<n; i++)
421 p[i] = pprsc(L); 421 p[i] = pprsc(L);
422 } 422 }
423 else 423 else
424 { 424 {
425 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */ 425 /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
426 const double sqrtL = sqrt(L); 426 const double sqrtL = sqrt(L);
427 for (i = 0; i < n; i++) 427 for (i = 0; i < n; i++)
428 { 428 {
429 p[i] = floor(RNOR*sqrtL + L + 0.5); 429 p[i] = floor(RNOR*sqrtL + L + 0.5);
430 if (p[i] < 0.0) 430 if (p[i] < 0.0)
431 p[i] = 0.0; /* will probably never happen */ 431 p[i] = 0.0; /* will probably never happen */
432 } 432 }
433 } 433 }
434 } 434 }
435 435
436 /* Generate one poisson variate */ 436 /* Generate one poisson variate */
437 double 437 double