Bug in rnorm. (PR#1664)

[Peter Dalgaard BSA]

Duncan Murdoch <dmurdoch@pair.com> writes:

> On Fri, 14 Jun 2002 00:39:58 +0200 (MET DST), p.dalgaard@biostat.ku.dk
> wrote:
> 
> >With your script, I get a nice horizontal line both with 1.4.0 and
> >1.5.0-patched.
> >
> >(RedHat Linux 7.1 on a PC) 
> 
> I get what Rolf described in 1.5.0-patched for Windows with the
> default generators: 
> 
> RNGkind()
> [1] "Marsaglia-Multicarry" "Kinderman-Ramage"  
> 
> I switched to the Ahrens-Dieter normal RNG (using
> RNGkind(,'Ahrens-Dieter')), and things were fine.  
> 
> Keeping the default normal generator but switching to the
> Mersenne-Twister uniform also fixed it:
> 
> > RNGkind()
> [1] "Mersenne-Twister" "Kinderman-Ramage"
> 
> Is it possible you're not using the default generators?  If you are,
> this is quite weird, because we're on the same hardware.  Why would
> the OS or compiler matter in a calculation like this??
> 
> Duncan Murdoch
> 

Now this is pretty darn odd....

This is what happened to me yesterday with 1.4.0:

[Previously saved workspace restored]

> m <- sapply(1:1000, function(i)max(rnorm(1500)^2))
> sum( 1 - pchisq(m,1)^1500 < .05)
[1] 39
> m <- sapply(1:1000, function(i)max(rnorm(1500)^2))
> sum( 1 - pchisq(m,1)^1500 < .05)
[1] 37
> m <- sapply(1:1000, function(i)max(rnorm(1500)^2))
> sum( 1 - pchisq(m,1)^1500 < .05)
[1] 63
> m <- sapply(1:1000, function(i)max(rnorm(1500)^2))
> sum( 1 - pchisq(m,1)^1500 < .05)
[1] 50
> m <- sapply(1:1000, function(i)max(rnorm(1500)^2))
> sum( 1 - pchisq(m,1)^1500 < .05)
[1] 51
> m <- sapply(1:1000, function(i)max(rnorm(1500)^2))
> sum( 1 - pchisq(m,1)^1500 < .05)
[1] 51

This morning, I get

> RNGkind()
[1] "Marsaglia-Multicarry" "Kinderman-Ramage"    
> m <- sapply(1:1000, function(i)max(rnorm(1500)^2))
> sum( 1 - pchisq(m,1)^1500 < .05)
[1] 27
> m <- sapply(1:1000, function(i)max(rnorm(1500)^2))
> sum( 1 - pchisq(m,1)^1500 < .05)
[1] 23
> m <- sapply(1:1000, function(i)max(rnorm(1500)^2))
> sum( 1 - pchisq(m,1)^1500 < .05)
[1] 29
> m <- sapply(1:1000, function(i)max(rnorm(1500)^2))
> sum( 1 - pchisq(m,1)^1500 < .05)
[1] 29

with the SAME version of R. The only difference is the startup
directory, and hence the saved workspace. This contains:

> ls(all=T)
[1] ".Random.seed" "x"            "y"           
> .Random.seed
[1]     0 25300 11635 10783

..whereas 1.5.0 which now again shows trouble has

>  .Random.seed
[1]         1 196857153 319227924

and 1.4.0 in a clean dir has this after a few iterations of the above

>  .Random.seed
[1]          1 1297109091   47580530


Oho! I think I get it: the "1" in the seed actually *defines* the RNG,
no matter what RNGkind() is telling me??

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