prev up next   top/contents search

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


prev up next   contents search
about this FAQ list   about eskimo   search   feedback   copyright

Hosted by Eskimo North