RandomLib::ExactNormal< bits > Class Template Reference

Sample exactly from a normal distribution. More...

#include <RandomLib/ExactNormal.hpp>

Public Member Functions
template<class Random >
RandomNumber< bits >	operator() (Random &r) const

Detailed Description

template<int bits = 1>
class RandomLib::ExactNormal< bits >

Sample exactly from a normal distribution.

Sample x from exp(−x²/2) / sqrt(2π). For background, see:

• J. von Neumann, Various Techniques used in Connection with Random Digits, J. Res. Nat. Bur. Stand., Appl. Math. Ser. 12, 36–38 (1951), reprinted in Collected Works, Vol. 5, 768–770 (Pergammon, 1963).
• D. E. Knuth and A. C. Yao, The Complexity of Nonuniform Random Number Generation, in "Algorithms and Complexity" (Academic Press, 1976), pp. 357–428.
• P. Flajolet and N. Saheb, The Complexity of Generating an Exponentially Distributed Variate, J. Algorithms 7, 463–488 (1986).

The algorithm is given in

• C. F. F. Karney, Sampling exactly from the normal distribution, http://arxiv.org/abs/1303.6257 (Mar. 2013).

In brief, the algorithm is:

1. Select an integer k ≥ 0 with probability exp(−k/2) (1−exp(−1/2)).
2. Accept with probability exp(− k (k − 1) / 2); otherwise, reject and start over at step 1.
3. Sample a random number x uniformly from [0,1).
4. Accept with probability exp(− x (x + 2k) / 2); otherwise, reject and start over at step 1.
5. Set x = k + x.
6. With probability 1/2, negate x.
7. Return x.

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:

• Step 1:
  • k = 0;
  • while (ExpProb(−1/2)) increment k by 1.
• Step 2:
  • n = k (k − 1) / 2;
  • while (n > 0) { if (!ExpProb(−1/2)) go to step 1; decrement n by 1; }
• Step 4:
  • repeat k + 1 times: if (!ExpProb(− x (x + 2k) / (2k + 2))) go to step 1.

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:

• Generate a sequence of random numbers U_i and find the greatest n such that 1/2 > U₁ > U₂ > . . . > U_n. (The resulting value of n may be 0.)
• If n is even, accept and return true; otherwise (n odd), reject and return false.

For p = x (x + 2k) / (2k + 2) (step 4), we generalize von Neumann's procedure as follows:

• Generate two sequences of random numbers U_i and V_i and find the greatest n such that both the following conditions hold
  • x > U₁ > U₂ > . . . > U_n;
  • V_i < (x + 2 k) / (2 k + 2) for all i in [1, n].
  (The resulting value of n may be 0.)
• If n is even, accept (return true); otherwise (n odd), reject (return false).

Here, instead of testing V_i < (x + 2 k) / (2 k + 2), we carry out the following tests:

• return true, with probability 2 k / (2 k + 2);
• return false, with probability 1 / (2 k + 2);
• otherwise (also with probability 1 / (2 k + 2)), return x > V_i.

The resulting method now entails evaluation of simple fractional probabilities (e.g., 1 / (2 k + 2)), or comparing random numbers (e.g., U₁ > U₂). 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.

#include <RandomLib/RandomNumber.hpp>
#include <RandomLib/ExactNormal.hpp>

  RandomLib::Random r;
  const int bits = 1;
  RandomLib::ExactNormal<bits> ndist;
  for (size_t i = 0; i < 10; ++i) {
    RandomLib::RandomNumber<bits> x = ndist(r); // Sample
    std::pair<double, double> z = x.Range();
    std::cout << x << " = "     // Print in binary with ellipsis
              << "(" << z.first << "," << z.second << ")"; // Print range
    double v = x.Value<double>(r); // Round exactly to nearest double
    std::cout << " = " << v << "\n";
  }

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.

Template Parameters

bits	the number of bits in each digit.

Definition at line 162 of file ExactNormal.hpp.

Member Function Documentation

template<int bits>

template<class Random >

RandomNumber< bits > RandomLib::ExactNormal< bits >::operator()

(

Random &

r

)

const

Return a random deviate with a normal distribution of mean 0 and variance 1.

Template Parameters

Random

the type of the random generator.

Parameters

[in,out]

r

a random generator.

Returns

the random sample.

Definition at line 308 of file ExactNormal.hpp.

---

The documentation for this class was generated from the following file:

• ExactNormal.hpp