3 Replies - 2071 Views - Last Post: 03 February 2012 - 05:10 PM

#1 jjhaag   User is offline

  • me editor am smartastic
  • member icon

Reputation: 48
  • View blog
  • Posts: 1,789
  • Joined: 18-September 07

Normally (Gaussian) Distributed Random Numbers

Posted 14 November 2007 - 12:35 PM

Description: Two function are provided, one that uses trig calls and one that does not - the latter is generally quicker. See the comments for a more detailed descriptionFunctions for creating normally-distributed random deviates, using the Box-Muller transformation.
/******************************************************************************/
/* randn()
 * 
 * Normally (Gaussian) distributed random numbers, using the Box-Muller 
 * transformation.  This transformation takes two uniformly distributed deviates
 * within the unit circle, and transforms them into two independently 
 * distributed normal deviates.  Utilizes the internal rand() function; this can
 * easily be changed to use a better and faster RNG.
 * 
 * The parameters passed to the function are the mean and standard deviation of 
 * the desired distribution.  The default values used, when no arguments are 
 * passed, are 0 and 1 - the standard normal distribution.
 * 
 * 
 * Two functions are provided:
 * 
 * The first uses the so-called polar version of the B-M transformation, using
 * multiple calls to a uniform RNG to ensure the initial deviates are within the
 * unit circle.  This avoids making any costly trigonometric function calls.
 * 
 * The second makes only a single set of calls to the RNG, and calculates a 
 * position within the unit circle with two trigonometric function calls.
 * 
 * The polar version is generally superior in terms of speed; however, on some
 * systems, the optimization of the math libraries may result in better 
 * performance of the second.  Try it out on the target system to see which
 * works best for you.  On my test machine (Athlon 3800+), the non-trig version
 * runs at about 3x10^6 calls/s; while the trig version runs at about
 * 1.8x10^6 calls/s (-O2 optimization).
 * 
 * 
 * Example calls:
 * randn_notrig();	//returns normal deviate with mean=0.0, std. deviation=1.0
 * randn_notrig(5.2,3.0);	//returns deviate with mean=5.2, std. deviation=3.0
 * 
 * 
 * Dependencies - requires <cmath> for the sqrt(), sin(), and cos() calls, and a
 * #defined value for PI.
 */

/******************************************************************************/
//	"Polar" version without trigonometric calls
double randn_notrig(double mu=0.0, double sigma=1.0) {
	static bool deviateAvailable=false;	//	flag
	static float storedDeviate;			//	deviate from previous calculation
	double polar, rsquared, var1, var2;
	
	//	If no deviate has been stored, the polar Box-Muller transformation is 
	//	performed, producing two independent normally-distributed random
	//	deviates.  One is stored for the next round, and one is returned.
	if (!deviateAvailable) {
		
		//	choose pairs of uniformly distributed deviates, discarding those 
		//	that don't fall within the unit circle
		do {
			var1=2.0*( double(rand())/double(RAND_MAX) ) - 1.0;
			var2=2.0*( double(rand())/double(RAND_MAX) ) - 1.0;
			rsquared=var1*var1+var2*var2;
		} while ( rsquared>=1.0 || rsquared == 0.0);
		
		//	calculate polar tranformation for each deviate
		polar=sqrt(-2.0*log(rsquared)/rsquared);
		
		//	store first deviate and set flag
		storedDeviate=var1*polar;
		deviateAvailable=true;
		
		//	return second deviate
		return var2*polar*sigma + mu;
	}
	
	//	If a deviate is available from a previous call to this function, it is
	//	returned, and the flag is set to false.
	else {
		deviateAvailable=false;
		return storedDeviate*sigma + mu;
	}
}


/******************************************************************************/
//	Standard version with trigonometric calls
#define PI 3.14159265358979323846

double randn_trig(double mu=0.0, double sigma=1.0) {
	static bool deviateAvailable=false;	//	flag
	static float storedDeviate;			//	deviate from previous calculation
	double dist, angle;
	
	//	If no deviate has been stored, the standard Box-Muller transformation is 
	//	performed, producing two independent normally-distributed random
	//	deviates.  One is stored for the next round, and one is returned.
	if (!deviateAvailable) {
		
		//	choose a pair of uniformly distributed deviates, one for the
		//	distance and one for the angle, and perform transformations
		dist=sqrt( -2.0 * log(double(rand()) / double(RAND_MAX)) );
		angle=2.0 * PI * (double(rand()) / double(RAND_MAX));
		
		//	calculate and store first deviate and set flag
		storedDeviate=dist*cos(angle);
		deviateAvailable=true;
		
		//	calcaulate return second deviate
		return dist * sin(angle) * sigma + mu;
	}
	
	//	If a deviate is available from a previous call to this function, it is
	//	returned, and the flag is set to false.
	else {
		deviateAvailable=false;
		return storedDeviate*sigma + mu;
	}
}


Is This A Good Question/Topic? 0
  • +

Replies To: Normally (Gaussian) Distributed Random Numbers

#2 Karel-Lodewijk   User is offline

  • D.I.C Addict
  • member icon

Reputation: 454
  • View blog
  • Posts: 864
  • Joined: 17-March 11

Re: Normally (Gaussian) Distributed Random Numbers

Posted 03 February 2012 - 05:08 PM

C++11 built-in support for random numbers with different distribution with among others normal distributed random numbers.
Was This Post Helpful? 0
  • +
  • -

#3 Karel-Lodewijk   User is offline

  • D.I.C Addict
  • member icon

Reputation: 454
  • View blog
  • Posts: 864
  • Joined: 17-March 11

Re: Normally (Gaussian) Distributed Random Numbers

Posted 03 February 2012 - 05:09 PM

The equivalent in c++11 would be std::normal_distribution distribution; //normal(0,1) std::mt19937 engine; // Mersenne twister MT19937 auto randn = std::bind(distribution, engine); engine.seed(time(0)); Then use randn() to get a normal(0,1) distributed number
Was This Post Helpful? 0
  • +
  • -

#4 Karel-Lodewijk   User is offline

  • D.I.C Addict
  • member icon

Reputation: 454
  • View blog
  • Posts: 864
  • Joined: 17-March 11

Re: Normally (Gaussian) Distributed Random Numbers

Posted 03 February 2012 - 05:10 PM

std::normal_distribution takes a template paramater float or double depending on the type of the random number to generate. The comment system just doesn't like the template brackets.
Was This Post Helpful? 0
  • +
  • -

Page 1 of 1