Actually, I loose few hours to search the code for random number generator which follows the Gaussian distribution with given mean and variance. Hence I am posting my final google and some hard work's result because I think this will help you to save more hours.
C++ Code (I have used Microsoft Visual Studio 2005). PLEASE PAY ATTENTION IN "stdafx.h" IF YOU ARE NOT USING Microsoft Visual Studio 2005:
// GussianDistributionGeneratior.cpp : Defines the entry point for the console application.
// Author: Rudra Poudel
#include "stdafx.h"
#include<stdio.h>
#include<stdlib.h>
#include<conio.h>
#include<math.h>
#include<float.h>
#include<limits.h>
#include<time.h>
//Random generator parameters
#define IA 16807
#define IM 2147483647
#define AM (1. /IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
long idum;
float getRand(long *idum);
float getGaussianRand(long *idum, double mean, double sigma);
int _tmain(int argc, _TCHAR* argv[])
{
idum= -time(0);
printf("Random Numbers:\n");
for(int i=0;i<100;i++){
printf("%f\t",getRand(&idum));
}
printf("\n\nPress any key to continue...");
getche();
printf("\n\nGaussian Random Numbers with mean=0 and variance=1:\n");
for(int i=0;i<100;i++){
printf("%f\t",getGaussianRand(&idum,4,1));
}
printf("\n\nPress any key to continue...");
getche();
return 0;
}
//Below function follows the Random Number Generators from Numerical Recipes in C, v2.0 (C) Copr. 1986-92 Numerical Recipes Software
/************************************************************************
* name: getRand & getGussianRand
* description: random generators from
* input: idum: Seed for random number generator
* output: Random number between 0-1 range
************************************************************************/
float getRand(long *idum)
{
int j;
long k;
static long iy=0;
static long iv[NTAB];
float temp;
if (*idum <= 0 || !iy) {
if (-(*idum) < 1) *idum = 1;
else *idum = -(*idum);
for (j=NTAB+7; j>=0; j--) {
k = (*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum <0) *idum += IM;
if (j < NTAB) iv[j] = *idum;
}
iy = iv[0];
}
k = (*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
j = iy/NDIV;
iy=iv[j];
iv[j] = *idum;
if ((temp=AM*iy) > RNMX) return RNMX;
else return temp;
}
/************************************************************************
* name: getGussianRand
* description: Gaussian random generators with given random number generator's seed, mean and variance
* input: idum: Seed for random number generator
* mean: mean for the distribution; for standard normal distribution mean = 0
* sigma: variance (square of standard deviation) for standard normal distribution sigma=1. As must of the case people need Gaussina Random number with mean =0 and sigma=1, hence if you are not sure just use that or plot the points and observer with different values.
* output: A random value form gussian distribution with given parameters
************************************************************************/
float getGaussianRand(long *idum, double mean, double sigma)
{
float getRand(long *idum);
static int iset=0;
static float gset;
float fac,rsq,v1,v2;
if (iset == 0) {
do {
v1=2.0*getRand(idum)-1.0;
v2=2.0*getRand(idum)-1.0;
rsq=v1*v1+v2*v2;
} while (rsq >= 1.0 || rsq==0.0);
fac=sqrt(-2.0*log(rsq)/rsq);
gset=v1*fac;
iset=1;
return ( (v2*fac*sigma) + mean);
} else {
iset =0;
return ( (gset*sigma) + mean);
}
}
As I am running out of time so I am unable to write description and theory behind it. Please don't hesitate to drop comment or question if you find something strange!!!
No comments:
Post a Comment