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