Generating random numbers from Normal distribution in C

Generating random numbers that are normally distributed

I needed to write a random number generator in C which will generate random numbers from Normal Distribution (Gaussian Distribution). Without this component I couldn’t proceed to finish writing a C code for Heuristic Kalman Algorithm by Lyonnet and Toscano for some experiments. I selected the Marsaglia and Bray method also known as the Polar method to generate Normal random variables. Here is how it is done.

I am just writing the algorithm

Generate $U1\thicksim\mathcal{U}(-1,1)$ and $U2\thicksim\mathcal{U}(-1,1)$
Generate $W$ until $W = U1^2 + U2^2 < 1$
Generate $X1=U1\cdot\sqrt{\frac{-2\ln(W)}{W}}$ and $X2=U2\cdot\sqrt{\frac{-2\ln(W)}{W}}$

Here $\mathcal{U}(-1,1)$ is the Uniform Distribution with range [-1,+1]
Now $X1$ and $X2$ are normal random variables with mean 0 and standard deviation 1. To generate normal random variable from mean $\mu$ and standard deviation $\sigma$ we need to do the following simple transform.

$X_{\mu,\sigma}=\mu+\sigma X_{0,1}$

Where $X_{0,1}\thicksim\mathcal{N}(0,1)$ and $X_{\mu,\sigma}\thicksim\mathcal{N}(\mu,\sigma)$

In each iteration two normal random variables are generated. Therefore we can generate two random variables in one iteration send one, and on the next call we will execute the algorithm and instead we will return the second generated value from the previous call. The implementation is pretty easy, the only thing we need is a uniform random number generator within the range [-1,+1]. We can use the uniform random number generator available in stblib, the rand function.

Sourcecode

Here is the code

#include <math.h>
#include <stdlib.h>

double
randn (double mu, double sigma)
{
  double U1, U2, W, mult;
  static double X1, X2;
  static int call = 0;

  if (call == 1)
    {
      call = !call;
      return (mu + sigma * (double) X2);
    }

  do
    {
      U1 = -1 + ((double) rand () / RAND_MAX) * 2;
      U2 = -1 + ((double) rand () / RAND_MAX) * 2;
      W = pow (U1, 2) + pow (U2, 2);
    }
  while (W >= 1 || W == 0);

  mult = sqrt ((-2 * log (W)) / W);
  X1 = U1 * mult;
  X2 = U2 * mult;

  call = !call;

  return (mu + sigma * (double) X1);
}

The rand () call returns a random number uniformly distributed within 0 to RAND_MAX. To generate uniform random numbers within range [0,1] we just need to divide the returned number with RAND_MAX. Here we need to explicitly typecast any of the operand to double to make the division floating point. Not doing it will result in an integer division which will always evaluate to 0 (and 1 if returned value is RAND_MAX). We need the uniform random number to be in range [-1,+1].

To scale a number in a range [low,high] we need the following transform, where x is scaled to range [low,high] and the scaled value is y

y = -low + x * (high - low)

Using the above transformation the statement -1 + ((double) rand () / RAND_MAX) * 2; generates uniform distribution in range [-1,+1]. X1 is normally distributed with mean 0 and standard deviation 1. We thus make the necessary transformation (mu + sigma * (double) X1) before returning the random variable, as in $X_{\mu,\sigma}=\mu+\sigma X_{0,1}$ .

The check for if W is 0, while (W >= 1 || W == 0) in the loop is done to avoid division by zero. The loop will keep on generating U1 and U2 as in the algorithm.

The variables X1 and X2 are made static so that it can hold the values from the previous call. Note that the value of X1 is not required to be held across iteration, but still is defined as static. The flag call determines if the call to the function randn is even or odd. If call = 0 then we generate two random numbers from normal distribution with mean 0 and standard deviation 1 using the Polar method, and then transform the generated random variable to make it have a mean mu and standard deviation sigma then return X1. If call = 1 then we do not compute anything and return the second normal random number (after mu sigma transformation) X2 generated in the previous call.