Bug in rnorm. (PR#1664)

Duncan Murdoch · 2002-06-13T22:39:58Z

Peter Dalgaard BSA writes: > rolf@math.unb.ca writes: > > > There appears to be a mild bug, or at least a deficiency, in > > rnorm. The bug becomes apparent when one looks at extremes > > of the squares of the values generated by rnorm; rnorm is not > > generating quite enough extreme values. > > > > The R version that I am using is 1.4.1; I never got around to installing > > 1.5.0, and now since 1.5.1 is about to come out .... However, checking > > the 1.5.0 rele

4 messages · Duncan Murdoch, Peter Dalgaard, Brian Ripley

Original

Peter Dalgaard

Thu, Jun 13, 2002 3:39 PM #

Peter Dalgaard BSA <p.dalgaard@biostat.ku.dk> writes:

No I do not! My quick check code took the square of the max rather
than the max of the square...

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)

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

Duncan Murdoch

Thu, Jun 13, 2002 6:03 PM #

On Fri, 14 Jun 2002 00:39:58 +0200 (MET DST), p.dalgaard@biostat.ku.dk
wrote:

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:

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

Peter Dalgaard

Fri, Jun 14, 2002 1:14 AM #

Duncan Murdoch <dmurdoch@pair.com> writes:

Now this is pretty darn odd....

This is what happened to me yesterday with 1.4.0:

[Previously saved workspace restored]

[1] 39

[1] 37

[1] 63

[1] 50

[1] 51

[1] 51

This morning, I get

[1] "Marsaglia-Multicarry" "Kinderman-Ramage"

[1] 27

[1] 23

[1] 29

[1] 29

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

[1] ".Random.seed" "x"            "y"

[1]     0 25300 11635 10783

..whereas 1.5.0 which now again shows trouble has

[1]         1 196857153 319227924

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

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

Brian Ripley

Fri, Jun 14, 2002 2:03 AM #

Just so we understand:

[1] "Marsaglia-Multicarry" "Kinderman-Ramage"

[1] "Marsaglia-Multicarry" "Kinderman-Ramage"

[1] 0.03381877

[1] "Wichmann-Hill"    "Kinderman-Ramage"

?RNGkind says

`RNGkind' returns a two-element character vector of the RNG and
normal kinds in use before the call,
                    ^^^^^^
However, any use of the random number generators themselves loads
.Random.seed (if it exist), and that will set the types of random
numbers in use.  So if a user changes .Random.seed (for example by
loading a saved workspace) the types will change.  RNGkind() reports
retrospectively (as it is documented to), not prospectively.

Now we could alter RNGkind() to read in the current .Random.seed (if there
is one) although it is a bit tricky to decide what to do if the value is
invalid.  (The random-number generation code throws an error and
expects thse user to fix, e.g  remove, .Random.seed.)

On 14 Jun 2002, Peter Dalgaard BSA wrote:

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??

Brian D. Ripley,                  ripley@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

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