I needed to write a random number generator in C which will generate random numbers from Normal Distribution (Gaussian Distribution). Without this component I couldn’t proceed to finish writing a C code for Heuristic Kalman Algorithm by Lyonnet and Toscano for some experiments. I selected the Marsaglia and Bray method also known as the Polar method to generate Normal random variables. Here is how it is done.

I am just writing the algorithm

- Generate and
- Generate until
- Generate and

Here is the Uniform Distribution with range `[-1,+1]`

Now and are normal random variables with mean 0 and standard deviation 1. To generate normal random variable from mean and standard deviation we need to do the following simple transform.

Where and

In each iteration two normal random variables are generated. Therefore we can generate two random variables in one iteration send one, and on the next call we will execute the algorithm and instead we will return the second generated value from the previous call. The implementation is pretty easy, the only thing we need is a uniform random number generator within the range `[-1,+1]`. We can use the uniform random number generator available in `stblib`, the `rand` function.

### Sourcecode

Here is the code

#include <math.h> #include <stdlib.h> double randn (double mu, double sigma) { double U1, U2, W, mult; static double X1, X2; static int call = 0; if (call == 1) { call = !call; return (mu + sigma * (double) X2); } do { U1 = -1 + ((double) rand () / RAND_MAX) * 2; U2 = -1 + ((double) rand () / RAND_MAX) * 2; W = pow (U1, 2) + pow (U2, 2); } while (W >= 1 || W == 0); mult = sqrt ((-2 * log (W)) / W); X1 = U1 * mult; X2 = U2 * mult; call = !call; return (mu + sigma * (double) X1); }

The `rand ()` call returns a random number uniformly distributed within `0` to `RAND_MAX`. To generate uniform random numbers within range `[0,1]` we just need to divide the returned number with `RAND_MAX`. Here we need to explicitly typecast any of the operand to `double` to make the division floating point. Not doing it will result in an integer division which will always evaluate to 0 (and 1 if returned value is `RAND_MAX`). We need the uniform random number to be in range `[-1,+1]`.

To scale a number in a range `[low,high]` we need the following transform, where `x` is scaled to range `[low,high]` and the scaled value is `y`

y = -low + x * (high - low)

Using the above transformation the statement `-1 + ((double) rand () / RAND_MAX) * 2;` generates uniform distribution in range `[-1,+1]`. `X1` is normally distributed with mean 0 and standard deviation 1. We thus make the necessary transformation `(mu + sigma * (double) X1)` before returning the random variable, as in .

The check for if `W` is 0, `while (W >= 1 || W == 0)` in the loop is done to avoid division by zero. The loop will keep on generating `U1` and `U2` as in the algorithm.

The variables `X1` and `X2` are made static so that it can hold the values from the previous call. Note that the value of `X1` is not required to be held across iteration, but still is defined as `static`. The flag `call` determines if the call to the function `randn` is even or odd. If `call = 0` then we generate two random numbers from normal distribution with mean 0 and standard deviation 1 using the Polar method, and then transform the generated random variable to make it have a mean `mu` and standard deviation `sigma` then return `X1`. If `call = 1` then we do not compute anything and return the second normal random number (after `mu` `sigma` transformation) `X2` generated in the previous call.

So that was all. Next I will post a simple fun code to plot histograms in terminal for this normal distribution (or any distribution or data set).

Generating random numbers with Gaussian (Normal) distribution in Visual Basic 6.0 (VB6) by Roberto Mior:

http://www.planetsourcecode.com/vb/scripts/ShowCode.asp?txtCodeId=75309&lngWId=1

Very useful for statistical analysis. Thank you very much indeed.

thanks a lot to Phoxis for help

please i need the simulation of normal distribution using monte carlo method can you help me

in C or C++

The values generated are always the same, how to solve? Sorry for my english!

Can you show me the code and how you are calling the function. I guess you are initialising the random number generator before each call.