Tuesday, March 27, 2012

Gaussian distribution random numbers generator in C

Hello!

Today a program to generate a sample of random numbers based on gaussian distribution.

To understand, you'll need to know C programming and basic random number theory.

Here is the code - explanations at the end.
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <stdint.h>
#include <sys/time.h>
#include <unistd.h>
#include <math.h>

main(int argc, char **argv)
{
struct timeval time;
int iRandom1, iRandom2;
double dRandom1, dRandom2, dAverage, dStddev;
long int i, liSampleSize;
   
//Incorrect number of parameters
if (argc!=4)
{
    printf ("Incorrect number of argument\n");
    printf ("Usage : %s SampleSize Average StdDev\n", argv[0]);
    printf ("Example : To have a 1000 samples centered with 1 of std deviation:\n");
    printf ("%s 1000 0 1\n", argv[0]);
    return 1;
}

//parsing the parameters
sscanf (argv[1],"%ld",&liSampleSize);
sscanf (argv[2],"%lf",&dAverage);
sscanf (argv[3],"%lf",&dStddev);


gettimeofday(&time, NULL);
   
//Convert seconds to a unsigned integer to init the seed
srandom((unsigned int) time.tv_usec);


for (i=0; i<liSampleSize;i++)
{
    //generates a number between 0 and 4G
    iRandom1=(random());
    //offset to have unifor ditribution between -2G and +2G
    dRandom1=(double)iRandom1 /2147483647;
   
    //generates a number between 0 and 4G
    iRandom2=(random());
    //offset to have unifor ditribution between -2G and +2G
    dRandom2=(double)iRandom2 /2147483647;
   
    //Output random values. The formula is based on Box-Muller method
    printf("%f \n", dAverage + dStddev * sqrt(-2.0 * log(dRandom1))*cos(6.28318531 * dRandom2));
}
return 0;
}


The code should be pasted in a file named gaussian.c
To compile the code:
gcc -lm gaussian.c  -o gaussian

To run it:
gaussian  500000 200 5  > dha.txt

This generates a list of 5000000 elements distributed around 200 with a std deviation of 5.

To check that you can either open dha.txt and check average and stdev in excel.

Or, on Linux you can run:
cat dha.txt | gsl-histogram 190 210 500 >dha_histogram.txt
awk '{print $1, $3 ; print $2, $3}' dha_histogram.txt | graph -T X
To display the distribution. You should get something similar to that
This is piece of code is fast efficient and will be useful for simulations in future.

Have fun!

1 comment: