From ripley at stats.ox.ac.uk  Tue May  2 18:54:20 2006
From: ripley at stats.ox.ac.uk (Prof Brian Ripley)
Date: Tue, 2 May 2006 17:54:20 +0100 (BST)
Subject: [Rd] normal random number generator in R
In-Reply-To: <200605021244.26160.tibshirani@gmail.com>
References: <200605021244.26160.tibshirani@gmail.com>
Message-ID: <Pine.LNX.4.64.0605021745570.24527@gannet.stats.ox.ac.uk>

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



[Rd] normal random number generator in R
Prof Brian Ripley