The discrete normal distribution. More...

#include <RandomLib/DiscreteNormal.hpp>

Public Member Functions
	DiscreteNormal (IntType sigma_num, IntType sigma_den=1, IntType mu_num=0, IntType mu_den=1)

template<class Random >
IntType	operator() (Random &r) const

Detailed Description

template<typename IntType = int>
class RandomLib::DiscreteNormal< IntType >

The discrete normal distribution.

Sample integers i with probability proportional to

\[ \exp\biggl[-\frac12\biggl(\frac{i-\mu}{\sigma}\biggr)^2\biggr], \]

where σ and μ are given as rationals (the ratio of two integers). The sampling is exact (provided that the random generator is ideal). For example

#include <iostream>
#include <RandomLib/Random.hpp>
#include <RandomLib/DiscreteNormal.hpp>
int main() {
  RandomLib::Random r;          // Create r
  r.Reseed();                   // and give it a unique seed
  int sigma_num = 7, sigma_den = 1, mu_num = 1, mu_den = 3;
  RandomLib::DiscreteNormal<int> d(sigma_num, sigma_den,
                                   mu_num, mu_den);
  for (int i = 0; i < 100; ++i)
    std::cout << d(r) << "\n";
}

prints out 100 samples with σ = 7 and μ = 1/3.

The algorithm is much the same as for ExactNormal; for details see

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

That algorithm samples the integer part of the result k, samples x in [0,1], and (unless rejected) returns s(k + x), where s = ±1. For the discrete case, we sample x in [0,1) such that

\[ s(k + x) = (i - \mu)/\sigma, \]

or

\[ x = s(i - \mu)/\sigma - k \]

The value of i which results in the smallest x ≥ 0 is

\[ i_0 = s\lceil k \sigma + s \mu\rceil \]

so sample

\[ i = i_0 + sj \]

where j is uniformly distributed in [0, ⌈σ⌉). The corresponding value of x is

\[ \begin{aligned} x &= \bigl(si_0 - (k\sigma + s\mu)\bigr)/\sigma + j/\sigma\\ &= x_0 + j/\sigma,\\ x_0 &= \bigl(\lceil k \sigma + s \mu\rceil - (k \sigma + s \mu)\bigr)/\sigma. \end{aligned} \]

After x is sampled in this way, it should be rejected if x ≥ 1 (this is counted with the next larger value of k) or if x = 0, k = 0, and s = −1 (to avoid double counting the origin). If x is accepted (in Step 4 of the ExactNormal algorithm), then return i.

When σ and μ are given as rationals, all the arithmetic outlined above can be carried out exactly. The basic rejection techniques used by ExactNormal are exact. Thus the result of this discrete form of the algorithm is also exact.

RandomLib provides two classes to sample from this distribution:

DiscreteNormal which is tuned for speed on a typical general purpose computer. This assumes that random samples can be generated relatively quickly.
DiscreteNormalAlt, which is a prototype for what might be needed on a small device used for cryptography which is using a hardware generator for obtaining truly random bits. This assumption here is that the random bits are relatively expensive to obtain.

The basic algorithm is the same in the two cases. The main advantages of this method are:

exact sampling (provided that the source of random numbers is ideal),
no need to cut off the tails of the distribution,
a short program involving simple integer operations only,
no dependence on external libraries (except to generate random bits),
no large tables of constants needed,
minimal time to set up for a new σ and μ (roughly comparable to the time it takes to generate one sample),
only about 5–20 times slower than standard routines to sample from a normal distribution using plain double-precision arithmetic.
DiscreteNormalAlt exhibits ideal scaling for the consumption of random bits, namely a constant + log₂σ, for large σ, where the constant is about 31.

The possible drawbacks of this method are:

σ and μ are restricted to rational numbers with sufficiently small numerators and denominators to avoid overflow (this is unlikely to be a severe restriction especially if the template parameter IntType is set to long long),
the running time is unbounded (but not in any practical sense),
the memory consumption is unbounded (but not in any practical sense),
the toll, about 30 bits, is considerably worse than that obtained using the Knuth-Yao algorithm, for which the toll is no more than 2 (but this requires a large table which is expensive to compute and requires a lot of memory to store).

This class uses a mutable private vector. So a single DiscreteNormal object cannot safely be used by multiple threads. In a multi-processing environment, each thread should use a thread-specific DiscreteNormal object.

Some timing results for IntType = int, μ = 0, and 10⁸ samples (time = time per sample, including setup time, rv = mean number of random variables per sample)

σ = 10, time = 219 ns, rv = 17.52
σ = 32, time = 223 ns, rv = 17.82
σ = 1000, time = 225 ns, rv = 17.95
σ = 160000, time = 226 ns, rv = 17.95

Template Parameters

IntType the integer type to use (default int).

Definition at line 140 of file DiscreteNormal.hpp.

Constructor & Destructor Documentation

template<typename IntType>

RandomLib::DiscreteNormal< IntType >::DiscreteNormal	(	IntType	sigma_num,
		IntType	sigma_den = `1`,
		IntType	mu_num = `0`,
		IntType	mu_den = `1`
	)

Constructor.

Parameters

[in]	sigma_num	the numerator of σ.
[in]	sigma_den	the denominator of σ (default 1).
[in]	mu_num	the numerator of μ (default 0).
[in]	mu_den	the denominator of μ (default 1).

The constructor creates a DiscreteNormal objects for sampling with specific values of σ and μ. This may throw an exception if the parameters are such that overflow is possible. Internally σ and μ are expressed with a common denominator, so it may be possible to avoid overflow by picking the fractions of these quantities so that sigma_den and mu_den have many common factors.

Definition at line 237 of file DiscreteNormal.hpp.

References STATIC_ASSERT.