[Rd] Replacing the Random Number Generator in Stand Alone Library

[peter dalgaard]

On Oct 16, 2013, at 02:41 , Nigel Delaney wrote:

> Okay, so I am guessing everyone had the same response I initially did when
> hearing that this RNG might not be so hot*...  as an alternate question, if
> I submitted a patch to replace the current RNG with the twister, would it be
> accepted?


Quite possibly. 

I think the reason you get no reply could be that nobody really knows to what extent the standalone library is being used, and what repercussions an internal change might have. (E.g., if a change in the seed format causes applications to crash and burn, users might not appreciate the improved RNG...)

The RNG issues really are serious, and affect actual applications. This was hashed out about 10 years ago and eventually fixed somewhere around R 1.6--1.7. 

> RNGversion("1.6.0")
Warning message:
In RNGkind("Marsaglia-Multicarry", "Buggy Kinderman-Ramage") :
  buggy version of Kinderman-Ramage generator used
> s <- replicate(1e6, max(rnorm(100)))
> plot(density(s))
> m <- matrix(runif(8e7),2)
> plot(m[1,],m[2,], xlim=c(0,1e-3), pch=".")

(The bug in the K-R generator isn't relevant here, but the pattern of RNG usage in K-R is what creates the ripple effect in the first plot. Changing _either_ of the two generators removes the effect.)

-pd


> 
> -N
> 
> -----Original Message-----
> From: Nigel Delaney [mailto:nigel.delaney at outlook.com] 
> Sent: Thursday, October 10, 2013 5:04 PM
> To: r-devel at r-project.org
> Subject: Replacing the Random Number Generator in Stand Alone Library
> 
> Hi R-Developers,
> 
> I had a question about the random number generator used in the R StandAlone
> Math Library.  The stand-alone library depends on the unif_rand() function
> for most simulated values, and this function is provided in the sunif.c file
> in the relevant directory.  At present, this program implements the
> "Marsaglia-Multicarry" algorithm, which is described throughout the R
> documentation as:
> 
> "A multiply-with-carry RNG is used, as recommended by George Marsaglia in
> his post to the mailing list 'sci.stat.math'. It has a period of more than
> 2^60 and has passed all tests (according to Marsaglia). The seed is two
> integers (all values allowed)."
> 
> However, I do not think this RNG actually passes all tests.   For example,
> the Handbook of Computational Econometrics (illegal web copy at link below),
> shows that it fails the mtuple test and gives an explicit example where it
> leads to problems because it failed this test.  The mtuple test was
> introduced by Marsaglia in 1985, and I gather he wrote his mailing list
> comment that it "passes all tests" sometime after this, so I am not sure
> what explains this distinction (though I am not sure if the mtuple test is
> included in the diehard tests, which he may have been what he was referring
> to).  However, there are clearly some areas where this PRNG runs in to
> trouble (although the books example is better, another problem is that it
> can't seem to simulate a value above (1/2)^1+(1/4)^4) after simulating a
> value below 1e-6.
> 
> Given that the Mersenne Twister seems to be the standard for simulation
> these days (and used as the default in R), it seems like it might be useful
> to change the stand alone library so it also uses this routine.  I gather
> this would be pretty easy to do by pulling this function from the RNG.c file
> and moving it into the sunif.c file, and have a prototype of this.
> 
> However, I thought I would ask, is there a reason this hasn't been done?  Or
> is it just a historical carry-over (pun intended I suppose).
> 
> Warm wishes,
> Nigel
> 
> ********************
> Research Fellow
> Massachusetts General Hospital / Broad Institute
> 
> 
> Book link:
> http://thenigerianprofessionalaccountant.files.wordpress.com/2013/04/handboo
> k-of-computational-econometrics-belsley.pdf
> 
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

-- 
Peter Dalgaard, Professor
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com