rnorm works correctly (was [R] rnorm??)

[Prof Brian Ripley]

On Mon, 21 Feb 2005, Scholz, Fritz wrote:

> I am wondering whether there is a bug in rnorm.

I am wondering if there is a bug in your procedures, or if you used a 
very old version of R.

If you call `bug', do have the courtesy to present your reasons.  See the 
posting guide for other missing information, like the R version, and in 
this case

> RNGkind()
[1] "Mersenne-Twister" "Inversion"

which shows the default.

> When generating rnorm(1000000) and counting
> the cases > 4 and the cases < (-4) I get rather
> unexpectedly low counts for the latter. The problem goes away
> when using qnorm(runif(1000000)).

That is the algorithm used by default internally (since 1.7.0, 
mid-2003, it seems).

> x <- rnorm(1000000); mean(x > 4); mean(x < -4)
[1] 3.5e-05
[1] 2.8e-05
> x <- rnorm(1000000); mean(x > 4); mean(x < -4)
[1] 3.5e-05
[1] 3.5e-05
> x <- rnorm(1000000); mean(x > 4); mean(x < -4)
[1] 2.7e-05
[1] 2.9e-05
> x <- rnorm(1000000); mean(x > 4); mean(x < -4)
[1] 4.3e-05
[1] 3.7e-05

looks about right.  As does

cnt <- 0
for(i in 1:100) {x <- rnorm(1000000); cnt <- cnt + sum(x < -4)}
cnt
[1] 3157

which is (to a good approximation) Poisson(3167)

> ppois(3157, 1e8*pnorm(-4))
[1] 0.4332376

or more exactly

> pbinom(3157, 1e8, pnorm(-4))
[1] 0.4332365


However, this is not a good way to generate such random events, and has 
many times tripped up people: see the cover of my 1987 simulation book for 
the theory behind one famous (and published) instance.

-- 
Brian D. Ripley,                  ripley at 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 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595