RandomLib
1.10

Sample exactly from a normal distribution. More...
#include <RandomLib/ExactNormal.hpp>
Public Member Functions  
template<class Random >  
RandomNumber< bits >  operator() (Random &r) const 
Sample exactly from a normal distribution.
Sample x from exp(−x^{2}/2) / sqrt(2π). For background, see:
The algorithm is given in
In brief, the algorithm is:
It is easy to show that this algorithm returns samples from the normal distribution with zero mean and unit variance. Futhermore, all these steps can be carried out exactly as follows:
Here, ExpProb(−p) returns true with probability exp(−p). With p = 1/2 (steps 1 and 2), this is implemented with von Neumann's rejection technique:
For p = x (x + 2k) / (2k + 2) (step 4), we generalize von Neumann's procedure as follows:
Here, instead of testing V_{i} < (x + 2 k) / (2 k + 2), we carry out the following tests:
The resulting method now entails evaluation of simple fractional probabilities (e.g., 1 / (2 k + 2)), or comparing random numbers (e.g., U_{1} > U_{2}). These may be carried out exactly with a finite mean running time.
With bits = 1, this consumes 30.1 digits on average and the result has 1.19 digits in the fraction. It takes about 676 ns to generate a result (1460 ns, including the time to round it to a double). With bits = 32, it takes 437 ns to generate a result (621 ns, including the time to round it to a double). In contrast, NormalDistribution takes about 44 ns to generate a double result.
Another way of assessing the efficiency of the algorithm is thru the mean value of the balance = (number of random bits consumed) − (number of bits in the result). If we code the result in Knuth & Yao's unarybinary notation, then the mean balance is 26.6.
For example the following samples from a normal exponential distribution and prints various representations of the result.
Here's a possible result:
1.00... = (1.25,1) = 1.02142 0.... = (1,0) = 0.319708 0.... = (0,1) = 0.618735 0.0... = (0.5,0) = 0.396591 0.0... = (0,0.5) = 0.20362 0.0... = (0,0.5) = 0.375662 1.111... = (2,1.875) = 1.88295 1.10... = (1.75,1.5) = 1.68088 0.... = (1,0) = 0.577547 0.... = (1,0) = 0.890553
First number is in binary with ... indicating an infinite sequence of random bits. Second number gives the corresponding interval. Third number is the result of filling in the missing bits and rounding exactly to the nearest representable double.
This class uses some mutable RandomNumber objects. So a single ExactNormal object cannot safely be used by multiple threads. In a multiprocessing environment, each thread should use a threadspecific ExactNormal object. In addition, these should be invoked with threadspecific random generator objects.
bits  the number of bits in each digit. 
Definition at line 162 of file ExactNormal.hpp.
RandomNumber< bits > RandomLib::ExactNormal< bits >::operator()  (  Random &  r  )  const 
Return a random deviate with a normal distribution of mean 0 and variance 1.
Random  the type of the random generator. 
[in,out]  r  a random generator. 
Definition at line 308 of file ExactNormal.hpp.