In order to produce normally distributed random numbers from the system
function rand() that outputs uniformly distributed random numbers, we will use
the "transformation method" described in W.H. Press, B.P. Flannery,
S.A. Teukolsky, and W.T. Vetterling, *Numerical Recipes in C: The Art of
Scientific Computing*, (Cambridge University Press, Cambridge, 1988). The
program segments below were adapted by Daniel Kiepfer from the
*Numerical Recipes* book.
The first listing is a "helper" function to produce a uniformly
distributed random number between two limits from the system supplied rand()
function.
//returns a uniformly distributed random number between min and max
double drandUniform (double min, double max)
{
return rand()*(max-min)/RAND_MAX+min;
}
Briefly, the idea of the transformation method is the following: If *f(x)*
is a uniform distribution, how is *f(y)* distributed given that *y=y(x)*?
The key is to examine the slope of the function *y(x)*, that is*,* to
look at *y’(x).* If *y’(x) = dy/dx* is large (in absolute value) for
some values of *x* and *y*, then *|dy| > |dx|.* The probability
of a number being in the range *dx* around *x* is the same for equal
length segments *dx*. That is, for a uniform distribution, a random number
has an equal probability of occurring *per unit x,* independent of *x.*
Upon transformation, the probability of occurring in the range *dx* about
*x* will be the same as the probability of occurring in the range *dy*
about *y*. Now if *|dy| > |dx|*, that means the probability *per
unit y* will be smaller. Also, of course, in general *y’(x)* will vary
with *x* and so *f(y)* will not be constant. One can show:
So, to produce a non-uniform distribution *f(y)*, one just needs to find
a function *y(x)* such that:
.
In principle this is all one needs, but the algorithm used below takes advantage
of a technique called the *Box-Muller* method, in which it turns out that
generating two random numbers with each function call takes less computational
effort than generating just one. This explains the structure of the function in
which alternate calls to the function return a value kept in a static variable
*gset*. Further details can be found in the *Numerical Recipes* book.
//returns a random number distributed normally, with a given SD
double drandGaussian (double stdDev)
{
static int noExtra = 1;
static double gset;
double fac, r, v1, v2;
if (noExtra)
{
do
{
v1 = drandUniform(-1.0, 1.0);
v2 = drandUniform(-1.0, 1.0);
r = v1*v1 + v2*v2;
} while (r >= 1.0);
fac = sqrt(-2.0*log(r)/r);
gset = v1*fac;
noExtra = 0;
return v2*fac*stdDev;
}
else
{
noExtra = 1;
return gset*stdDev;
}
} |