Normal Distribution Random Number

Hello,

A good method to generate random numbers with a probability distribution y = f(x) is as in the steps given:

1. Seed the random number. One can use the function srand(time(NULL)) for this. It is best to introduce in the function which generates random numbers, to be run at the first time.

2. One can write a function which generates random numbers between values x0 and x1 as follows:

#include <time.h>
#include <stdlib.h>

double Random(double x0, double x1)
{
static bool Flag = false;
if (!Flag)
{
Flag = true;
srand((unsigned) time(NULL));
}
return x0 + (x1 - x0) * rand() / RAND_MAX;
}

Note that the seed generation is embedded in the function implementation and will be called the first time.

3. Determine the maximum y value of the function within the range of interest. Let it be y1. Minimum y value will normally be y0 = 0.

4. Generate two random numbers x = Random(x0, x1) and y = Random(y0, y1). Check whether y < f(x). Accept the random number x if the condition is true, discard it if it is not.

5. The set of random numbers thus generated will be as per the distribution function y = f(x).

For probability functions which extends upto infinity on one or both sides of the limits, one has to map it to finite limits by changing variable. x=+/- (infinity) can be mapped to z=+/- (Pi/2) if we change the variable to x=cos(z).

Best luck.
Pravin.
10-03-2004.