Prof Darren Wilkinson Professor of Stochastic Modelling
School of Mathematics & Statistics
Newcastle University

Normal random variables via the Box-Muller method

We know that if T is Unif(0,2*Pi) and R^2 is Exp(0.5), then X=R*cos(T) is N(0,1). Since we know how to generate uniforms and exponentials, we can use this to generate standard Normal random quantities.

Implementation in R

We can write a vectorised R function to do this as follows. 

myrnorm<-function (n) 
{
        theta <- runif(n, 0, 2 * pi)
        rsq <- rexp(n, 0.5)
        x <- sqrt(rsq) * cos(theta)
        x
}
theta is a vector of n Unif(0,2*Pi) random quantities, rsq is a vector of Exp(0.5) random quantities, and so we use these to construct x to be our vector of N(0,1) random quantities.

The R source code for this is here: myrnorm.r

Again, because this procedure can be entirely vectorised, it will be very fast and efficient in R, and so there will be little to be gained from re-coding in C.

Implementation in C

In order to get a bit more of a feel for C, we'll see how to write some code to generate Normal random quantities. First rename the urand function from before to be surand (standard-urand). Now define a more general urand as follows. 
float urand(float low, float high)
{
  return(low+(high-low)*surand());
}
so now urand(a,b) will generate a Unif(a,b) random quantity. The genexp function from before will now need to be modified to call surand rather than urand. We can now define a standard Normal generation function as follows. 
float gennor()
{
  float theta,rsq,x;
  theta=urand(0,2*PI);
  rsq=genexp(0.5);
  x=sqrt(rsq)*cos(theta);
  return(x);
}

Full C source code is here: norm.c

Other implementations

Those who are interested may also wish to look at corresponding code for LISP-STAT: mynorm-rand.lsp, Python: norm.py and Sather: norm.sa.

Darren J Wilkinson Book: Stochastic Modelling for Systems Biology  
darren.wilkinson@ncl.ac.uk
http://www.staff.ncl.ac.uk/d.j.wilkinson/