R devel - May 2006 - normal random number generator in R

Hi,

I am trying to figure out how normal random number generator works in R. As I 
look at .../src/nmath/snorm.c file, I find the default algorithm is inverse 
CDF method. In more detail, instead of directly using uniform value by 
unif_rand(), snorm function will first get a sum by adding unif_rand() and 
2^27*unif_rand(), then divide it by 2^27 and transfer it to qnorm5() function 
for inverting. Only a short comment for this operation is available in the 
source:
/* unif_rand() alone is not of high enough precision */

Just curious why this operation is needed? Is it a general algorithm for 
inverse CDF method, or simply unif_rand() in R returns float precision?

Thanks
Tib

On Tue, 2 May 2006, Tibshirani wrote:
> Hi,
>
> I am trying to figure out how normal random number generator works in R. As
I
> look at .../src/nmath/snorm.c file, I find the default algorithm is inverse
> CDF method. In more detail, instead of directly using uniform value by
> unif_rand(), snorm function will first get a sum by adding unif_rand() and
> 2^27*unif_rand(), then divide it by 2^27 and transfer it to qnorm5()
function
> for inverting. Only a short comment for this operation is available in the
> source:
> /* unif_rand() alone is not of high enough precision */
>
> Just curious why this operation is needed? Is it a general algorithm for
> inverse CDF method, or simply unif_rand() in R returns float precision?
unif_rand() returns a uniform over about 2^31 distinct doubles, and that 
gives granularity in the far tails of a normal.

It is a known (at least to me) trick for improving the accuracy of 
inversion.

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