GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
randpoisson.c
Go to the documentation of this file.
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_H)
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 NAN
48 #define NAN octave_NaN
49 #undef INFINITE
50 #define INFINITE lo_ieee_isinf
51 #define RUNI oct_randu()
52 #define RNOR oct_randn()
53 #define LGAMMA xlgamma
54 
56 F77_FUNC (dlgams, DLGAMS) (const double *, double *, double *);
57 
58 static double
59 xlgamma (double x)
60 {
61  double result;
62 #ifdef HAVE_LGAMMA
63  result = lgamma (x);
64 #else
65  double sgngam;
66 
67  if (lo_ieee_isnan (x))
68  result = x;
69  else if (x <= 0 || lo_ieee_isinf (x))
70  result = octave_Inf;
71  else
72  F77_XFCN (dlgams, DLGAMS, (&x, &result, &sgngam));
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 C0 9.18938533204672742e-01
84 #define C1 8.33333333333333333e-02
85 #define C3 -2.77777777777777778e-03
86 #define C5 7.93650793650793651e-04
87 #define C7 -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 + C0 + r*(C1 + rr*(C3 + rr*(C5 + rr*C7))));
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 = RUNI * 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 * RUNI) + 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 * RUNI) + 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 = RUNI;
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 TABLESIZE 46
316  double t[TABLESIZE];
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 = RUNI;
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 < TABLESIZE)
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[TABLESIZE];
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 = RUNI;
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 < TABLESIZE)
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 - LGAMMA(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_PI*RUNI);
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 (RUNI > 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 - LGAMMA(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_PI*RUNI);
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 (RUNI > 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))
486  {
487  for (i=0; i<n; i++)
488  p[i] = 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 (RNOR*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 = 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 *= RUNI;
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))
537  {
538  /* FIXME R uses NaN, but the normal approx. suggests that as
539  * limit should be inf. Which is correct? */
540  ret = NAN;
541  }
542  else
543  {
544  /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
545  ret = floor (RNOR*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))
558  {
559  for (i=0; i<n; i++)
560  p[i] = 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 (RNOR*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 = 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 *= RUNI;
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))
610  {
611  /* FIXME R uses NaN, but the normal approx. suggests that as
612  * limit should be inf. Which is correct? */
613  ret = NAN;
614  }
615  else
616  {
617  /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
618  ret = floor (RNOR*sqrt (L) + L + 0.5);
619  if (ret < 0.0) ret = 0.0; /* will probably never happen */
620  }
621  return ret;
622 }