Strange means of numbers drawn from rpois (PR#730)

Kjetil Kjernsmo <kjetil.kjernsmo@astro.uio.no> writes:

On 8 Nov 2000, Peter Dalgaard BSA wrote:

range(sapply(1:2000, function(n) mean(rpois(10000, c(15,15+1e-8)))))

[1] 14.8692 15.1200

AHA! Spotted, I think.

Wow! Great, that was fast!

It is possible to return from rpois in step N,
in which case the initialisations in step P are not performed in that
invokation and not the next times either because it thinks that it has
been done previously.

Ah, I tried to poke at the code, but got confused rather fast...

I'm not quite sure what is the best way to fix it, though. Either kill
muprev when exiting in step N, or make the initialisations currently
in step P before entering step N. Or maybe save two versions of muprev
... yes, that might be it - will try.

Great! Thanks, I hope it'll work out. And thanks to all who responded!

Seems to have done the trick:

sapply(1:20,function(n) mean(rpois(100000, n)) )

[1]  1.00021  1.99588  2.99987  4.00364  4.98532  6.00528  7.00402  8.00046
 [9]  8.98958 10.00336 10.99647 12.00490 13.00697 14.01364 14.99357 15.99705
[17] 16.97901 17.99855 19.00363 20.02180


This is the patch (will commit to r-devel shortly)

[pd@blueberry R]$ cvs diff -u src/nmath/rpois.c 
Index: src/nmath/rpois.c
===================================================================
RCS file: /home/rdevel/CVS-ARCHIVE/R/src/nmath/rpois.c,v
retrieving revision 1.11
diff -u -r1.11 rpois.c
--- src/nmath/rpois.c   2000/08/30 16:45:16     1.11
+++ src/nmath/rpois.c   2000/11/08 18:09:34
@@ -65,7 +65,7 @@
 
     static double big_l;/* integer "w/o overflow" */
     static int l, m;
-    static double muprev = 0.;/*, muold         = 0.*/
+    static double muprev = 0., muprev2 = 0.;/*, muold   = 0.*/
 
     /* Local Vars  [initialize some for -Wall]: */
     double del, difmuk= 0., E= 0., fk= 0., fx, fy, g, px, py, t, u= 0., v, x;
@@ -160,7 +160,11 @@
     /* Step P. preparations for steps Q and H.
        (recalculations of parameters if necessary) */
 
-    if (new_big_mu) {
+    if (new_big_mu || mu != muprev2) {
+        /* Careful! muprev2 is not always == muprev
+          because one might have exited in step I or S
+          */
+        muprev2 = mu; 
        omega = M_1_SQRT_2PI / s;
        /* The quantities b1, b2, c3, c2, c1, c0 are for the Hermite
         * approximations to the discrete normal probabilities fk. */

O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk)             FAX: (+45) 35327907

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._