unexpected behaviour of rnorm(): summary

Thanks to everyone who contributed to the discussion on rnorm(), both
online and off.  I learned a lot.  Several people pointed out that the
issue was reported in bug report #1664 and indeed it was.  Also, the
problem goes away with RNGkind(, "Inversion").

However, I regard trying to break a random number generator as
precisely the best way to understand random variables and hypothesis
testing (does anyone else agree? Has anyone else tried to break
rnorm() with or without success? can people post their attempts?)

The original script was

f <- function(n){max(rnorm(n))}
plot(sapply(rep(5000,4000),f))   #[this takes my PC about 30 seconds]

Received wisdom is that "The Marsaglia-Multicarry and Kinderman-Ramage
options don't play nicely together in the extreme tails of the Normal
distribution" (PD) but I discovered last night that

f.uppertail <- function(n,upper=3){x <- rnorm(n);x[x>upper]}
plot(f.uppertail(1e7))

shows the effect much more directly...and viewed this way, the default
RNG would seem to be more seriously flawed.  The second script above
considers the upper (ie Z>3) tail of the Gaussian: and it has gaps.
This is borderline "extreme tails" because pnorm(3) is only about
0.998.

In contrast, Bug #1664 refers to the maximal value of a series of
rnorm() values which, as PBR points out, _is_ a rather severe test of
any RNG; the difference between f.uppertail() and #1664 would be that
#1664 simulates

dmaxnorm <- function(x,n){n* dnorm(x)*pnorm(x)^(n-1)}

which involves PDFs and CDFs; f.uppertail() is solely a rnorm() issue.
If nothing else, f.uppertail() is a simpler test for rnorm().

many thanks again to the List--it Rocks

rksh

Robin Hankin, Lecturer,
School of Geography and Environmental Science
Tamaki Campus
Private Bag 92019 Auckland
New Zealand

r.hankin at auckland.ac.nz
tel 0064-9-373-7599 x6820; FAX 0064-9-373-7042

as of: Fri Nov 29 09:27:00 NZDT 2002
This (linux) system up continuously for:  456 days, 15 hours, 09 minutes
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help 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-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._