[R] random number generator?

[ripley@stats.ox.ac.uk]

Can I suggest

RNGkind("Mersenne-Twister", "Inversion")

and especially the use of Inversion where tail behaviour of the normal is
important.

Were it not for concerns about reproducibility we would have switched to 
Inversion a while back.

On Tue, 28 Jan 2003, Charles Annis, P.E. wrote:

> 
> Earlier today I reported finding an unbalanced number of observations in
> the p=0.0001 tails of rnorm.
> 
> Many thanks to Peter Dalgaard who suggested changing the normal.kind
> generator.  
> 
> Using  RNGkind(kind = NULL, normal.kind ="Box-Muller")
> seems to have provided the remedy.  For example:
> 
> > observed.fraction.below 
> [1] 0.000103
> > observed.fraction.above 
> [1] 0.000101
> >  
> 
> Thank you, Peter!
> 
> 
> Charles Annis, P.E.
> 
> Charles.Annis at StatisticalEngineering.com
> phone: 561-352-9699
> eFAX: 503-217-5849
> http://www.StatisticalEngineering.com
> 
> 
> -----Original Message-----
> From: r-help-admin at stat.math.ethz.ch
> [mailto:r-help-admin at stat.math.ethz.ch] On Behalf Of Peter Dalgaard BSA
> Sent: Tuesday, January 28, 2003 2:36 PM
> To: Charles Annis, P.E.
> Cc: r-help at stat.math.ethz.ch
> Subject: Re: [R] random number generator?
> 
> "Charles Annis, P.E." <AnnisC at asme.org> writes:
> 
> > Dear R-Aficionados:
> > 
> > I realize that no random number generator is perfect, so what I report
> > below may be a result of that simple fact.  However, if I have made an
> > error in my thinking I would greatly appreciate being corrected.
> > 
> > I wish to illustrate the behavior of small samples (n=10) and so
> > generate 100,000 of them.
> > 
> > n.samples <- 1000000
> > sample.size = 10
> > p <- 0.0001
> > z.normal <- qnorm(p)
> > # generate n.samples of sample.size each from a normal(mean=0, sd=1)
> > density
> > #
> > small.sample <- matrix(rnorm(n=sample.size*n.samples, mean=0, sd=1),
> > nrow=n.samples, ncol=sample.size)
> > # Verify that from the entire small.sample matrix, p sampled values
> are
> > below, p above.
> > #
> > observed.fraction.below <- sum(small.sample <
> > z.normal)/length(small.sample)
> > observed.fraction.above <- sum(small.sample >
> > -z.normal)/length(small.sample)
> > 
> > > observed.fraction.below 
> > [1] 6.3e-05
> > > observed.fraction.above 
> > [1] 0.000142
> > > 
> > 
> > I've checked the behavior of the entire sample's mean and median and
> > they seem fine.  The total fraction in both tails is 0.0002, as it
> > should be.  However in every instance about 1/3 are in the lower tail,
> > 2/3 in the upper.  I also observe the same 1/3:2/3 ratio for one
> million
> > samples of ten.
> > 
> > Is this simply because random number generators aren't perfect?  Or
> have
> > I stepped in something?
> > 
> > Thank you for your kind counsel.
> 
> You stepped in something, I think, but I probably shouldn't elaborate
> on the metaphor ... There's an unfortunate interaction between the two
> methods that are used for generating uniform and normal variables (the
> latter uses the former). This has been reported a couple of times
> before and typically gives anomalous tail behaviour. Changing one of
> the generators (see help(RNGkind)) usually helps.
> 
> 

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