Box-Muller transform (C)

From LiteratePrograms

Suppose our random number generator is rand() and returns an integer between 0 and RAND_MAX, inclusive; any uniform random number generator will suffice. We produce two random integers and scale them to floating-point numbers between −1 and 1:

<<scale two random integers to doubles between -1 and 1>>=
x = 2.0*rand()/RAND_MAX - 1;
y = 2.0*rand()/RAND_MAX - 1;

To use the standard library's rand() we add stdlib.h:

<<header files>>=
#include <stdlib.h>

Diagram of the Box Muller transform. The initial circles, uniformly spaced about the origin, are mapped to another set of circles about the origin that are closely spaced near the origin but quickly spread out. The largest circles in the domain map to the smallest circles in the range and vice versa.

In order to use the polar form of the Box-Muller transform, it is necessary that the point (x,y) falls inside the unit circle (the circle of radius 1 about the origin), and not at the origin. To ensure this, we simply reject any pairs that don't meet these criteria and try again:

<<choose a point x,y in the unit circle uniformly at random>>=
double x, y, r;
do {
    scale two random integers to doubles between -1 and 1
    r = x*x + y*y;
} while (r == 0.0 || r > 1.0);

Above, r is just the square of the Euclidean distance of the point (x,y) from the origin (0,0). We don't need to take the square root because if the square falls between 0 and 1, so will the distance, and likewise if it is greater than 1 (this is a property of the square root function).

Next, we compute the desired normally distributed value pair using the Box Muller transform:

<<Apply Box-Muller transform on x, y>>=
double d = sqrt(-2.0*log(r)/r);
double n1 = x*d;
n2 = y*d;

This is the most complex part to understand. Intuitively, it maps each circle of points around the origin to another circle of points around the origin, as shown in the diagram, where larger outer circles are mapped to closely-spaced inner circles and inner circles to outer circles. The sqrt() and log() functions require a header:

<<header files>>=
#include <math.h>

These points will be distributed with a mean of 0 and a standard deviation of 1. To achieve any mean and standard deviation, we simply scale and translate:

<<scale and translate to get desired mean and standard deviation>>=
double result = n1*stddev + mean;

Some applications may want both values returned; however, since most will want only a single value at a time, we simply return one of the values and cache the other one for the next call, like this:

<<complete Box-Muller function>>=
double rand_normal(double mean, double stddev) {
    static double n2 = 0.0;
    static int n2_cached = 0;
    if (!n2_cached) {
        choose a point x,y in the unit circle uniformly at random
        {
        Apply Box-Muller transform on x, y
        scale and translate to get desired mean and standard deviation
        n2_cached = 1;
        return result;
        }
    } else {
        n2_cached = 0;
        return n2*stddev + mean;
    }
}

Alternatively we could just drop the second value, but the sqrt and log values are expensive enough that it's worth calling them less often.

We can now call the code with a simple test driver like this:

<<box_muller.c>>=
header files
#include <time.h>
#include <stdio.h>
complete Box-Muller function
int main() {
    int buckets[20] = {0};
    int i, j;
    srand(time(NULL));
    for(i=0; i<500; i++) {
        double val = rand_normal(10.0, 3.0);
        int i = (int)floor(val + 0.5);
        if (i >= 0 && i < sizeof(buckets)/sizeof(*buckets))
            buckets[i]++;
    }
    for(i=0; i<sizeof(buckets)/sizeof(*buckets); i++) {
        for(j=0; j<buckets[i]; j++) {
            printf("#");
        }
        printf("\n");
    }
    return 0;
}

The test driver creates a simple histogram of the output values by incrementing the "bucket" that each value falls into. After doing enough of these, the shape should appear roughly like a normal curve. Here's some sample output:

#
###
############
################
########################
##################################################
###########################################################
##############################################################
####################################################################
##########################################################
###############################################
############################################
########################
#################
########
######
#