[Rd] Re: [R] Strange means of numbers drawn from rpois (PR#730)

[p.dalgaard@biostat.ku.dk]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._