Nigel Delaney
2013-Oct-10 21:04 UTC
[Rd] 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
Nigel Delaney
2013-Oct-16 00:41 UTC
[Rd] FW: Replacing the Random Number Generator in Stand Alone Library
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? -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