Question 13.20

comp.lang.c FAQ list · Question 13.20

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

---

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.Additional links
References: Knuth Sec. 3.4.1 p. 117
Box and Muller, ``A Note on the Generation of Random Normal Deviates''
Marsaglia and Bray, ``A Convenient Method for Generating Normal Variables''
Abramowitz and Stegun, Handbook of Mathematical Functions
Press et al., Numerical Recipes in C Sec. 7.2 pp. 288-290