Sample exactly from a normal distribution. More...
|template<class Random >|
|RandomNumber< bits >||operator() (Random &r) const|
Sample exactly from a normal distribution.
Sample x from exp(−x2/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 Vi < (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., U1 > U2). 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 unary-binary 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 multi-processing environment, each thread should use a thread-specific ExactNormal object. In addition, these should be invoked with thread-specific random generator objects.
|bits||the number of bits in each digit.|
|RandomNumber< bits > RandomLib::ExactNormal< bits >::operator()||(||Random &||r||)||const|