How can I generate random numbers with a normal or Gaussian distribution?

Q

How can I generate random numbers with a normal or Gaussian distribution?

✍: Guest

A

There are a number of ways of doing this.
1. Exploit the Central Limit Theorem (``law of large numbers'') and add up several uniformly-distributed random numbers:
#include <stdlib.h>
#include <math.h>

#define NSUM 25

double gaussrand()
{
double x = 0;
int i;
for(i = 0; i < NSUM; i++)
x += (double)rand() / RAND_MAX;

x -= NSUM / 2.0;
x /= sqrt(NSUM / 12.0);

return x;
}

(Don't overlook the sqrt(NSUM / 12.) correction, though it's easy to do so accidentally, especially when NSUM is 12.)
2. Use a method described by Abramowitz and Stegun:

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

#define PI 3.141592654

double gaussrand()
{
static double U, V;
static int phase = 0;
double Z;

if(phase == 0) {
U = (rand() + 1.) / (RAND_MAX + 2.);
V = rand() / (RAND_MAX + 1.);
Z = sqrt(-2 * log(U)) * sin(2 * PI * V);
} else
Z = sqrt(-2 * log(U)) * cos(2 * PI * V);

phase = 1 - phase;

return Z;
}

3. Use a method discussed in Knuth and due originally to Marsaglia:

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

double gaussrand()
{
static double V1, V2, S;
static int phase = 0;
double X;

if(phase == 0) {
do {
double U1 = (double)rand() / RAND_MAX;
double U2 = (double)rand() / RAND_MAX;

V1 = 2 * U1 - 1;
V2 = 2 * U2 - 1;
S = V1 * V1 + V2 * V2;
} while(S >= 1 || S == 0);

X = V1 * sqrt(-2 * log(S) / S);
} else
X = V2 * sqrt(-2 * log(S) / S);

phase = 1 - phase;

return X;
}

These methods all generate numbers with mean 0 and standard deviation 1. (To adjust to some other distribution, multiply by the standard deviation and add the mean.) Method 1 is poor ``in the tails'' (especially if NSUM is small), but methods 2 and 3 perform quite well. See the references for more information.

2015-07-22, 903👍, 0💬