[R] random number generator?
Bob Wheeler
bwheeler at echip.com
Wed Jan 29 15:11:05 CET 2003
Well if the base generator is to be changed, it might be a
good idea to consider implementing one of Marsaglia's really
long period generators. His latest, using titanic primes,
has a period of >2^131086. It is extremely fast. Even more
could be done, by replacing the two phase normal generator
by faster generator using a single phase.
"Liaw, Andy" wrote:
>
> Might I suggest taking a poll (even though unscientific) of how many people
> will be affected by a change in default RNG? My totally arbitrary guess is
> very few, if any.
>
> If I'm not mistaken, Python had only recently changed the default RNG to
> Mersenne-Twister. If Python can do it, I should think R can, too, without
> too much pain...
>
> Just my $0.02...
>
> Andy
>
> > -----Original Message-----
> > From: ripley at stats.ox.ac.uk [mailto:ripley at stats.ox.ac.uk]
> > Sent: Tuesday, January 28, 2003 3:53 PM
> > To: Charles Annis, P.E.
> > Cc: r-help at stat.math.ethz.ch
> > Subject: RE: [R] random number generator?
> >
> >
> > 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
> >
> > ______________________________________________
> > R-help at stat.math.ethz.ch mailing list
> > http://www.stat.math.ethz.ch/mailman/listinfo/r-help
> >
>
> ------------------------------------------------------------------------------
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> http://www.stat.math.ethz.ch/mailman/listinfo/r-help
--
Bob Wheeler --- http://www.bobwheeler.com/
ECHIP, Inc. --- (302) 239-6620, voice FAX
724 Yorklyn Rd., Hockessin, DE 19707
Randomness comes in bunches
More information about the R-help
mailing list