RandomLib  1.10
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
NormalDistribution.hpp
Go to the documentation of this file.
1 /**
2  * \file NormalDistribution.hpp
3  * \brief Header for NormalDistribution
4  *
5  * Compute normal deviates.
6  *
7  * Copyright (c) Charles Karney (2006-2011) <charles@karney.com> and licensed
8  * under the MIT/X11 License. For more information, see
9  * http://randomlib.sourceforge.net/
10  **********************************************************************/
11 
12 #if !defined(RANDOMLIB_NORMALDISTRIBUTION_HPP)
13 #define RANDOMLIB_NORMALDISTRIBUTION_HPP 1
14 
15 #include <cmath> // for std::log
16 
17 namespace RandomLib {
18  /**
19  * \brief Normal deviates
20  *
21  * Sample from the normal distribution.
22  *
23  * This uses the ratio method; see Knuth, TAOCP, Vol 2, Sec. 3.4.1.C,
24  * Algorithm R. Unlike the Box-Muller method which generates two normal
25  * deviates at a time, this method generates just one. This means that this
26  * class has no state that needs to be saved when checkpointing a
27  * calculation. Original citation is\n A. J. Kinderman, J. F. Monahan,\n
28  * Computer Generation of Random Variables Using the Ratio of Uniform
29  * Deviates,\n ACM TOMS 3, 257--260 (1977).
30  *
31  * Improved "quadratic" bounds are given by\n J. L. Leva,\n A Fast Normal
32  * Random Number Generator,\n ACM TOMS 18, 449--453 and 454--455
33  * (1992).
34  *
35  * The log is evaluated 1.369 times per normal deviate with no bounds, 0.232
36  * times with Knuth's bounds, and 0.012 times with the quadratic bounds.
37  * Time is approx 0.3 us per deviate (1GHz machine, optimized, RealType =
38  * float).
39  *
40  * Example
41  * \code
42  * #include <RandomLib/NormalDistribution.hpp>
43  *
44  * RandomLib::Random r;
45  * std::cout << "Seed set to " << r.SeedString() << "\n";
46  * RandomLib::NormalDistribution<double> normdist;
47  * std::cout << "Select from normal distribution:";
48  * for (size_t i = 0; i < 10; ++i)
49  * std::cout << " " << normdist(r);
50  * std::cout << "\n";
51  * \endcode
52  *
53  * @tparam RealType the real type of the results (default double).
54  **********************************************************************/
55  template<typename RealType = double> class NormalDistribution {
56  public:
57  /**
58  * The type returned by NormalDistribution::operator()(Random&)
59  **********************************************************************/
60  typedef RealType result_type;
61  /**
62  * Return a sample of type RealType from the normal distribution with mean
63  * &mu; and standard deviation &sigma;.
64  *
65  * For &mu; = 0 and &sigma; = 1 (the defaults), the distribution is
66  * symmetric about zero and is nonzero. The maximum result is less than 2
67  * sqrt(log(2) \e p) where \e p is the precision of real type RealType.
68  * The minimum positive value is approximately 1/2<sup><i>p</i>+1</sup>.
69  * Here \e p is the precision of real type RealType.
70  *
71  * @tparam Random the type of RandomCanonical generator.
72  * @param[in,out] r the RandomCanonical generator.
73  * @param[in] mu the mean value of the normal distribution (default 0).
74  * @param[in] sigma the standard deviation of the normal distribution
75  * (default 1).
76  * @return the random sample.
77  **********************************************************************/
78  template<class Random>
79  RealType operator()(Random& r, RealType mu = RealType(0),
80  RealType sigma = RealType(1)) const throw();
81  };
82 
83  template<typename RealType> template<class Random> inline RealType
84  NormalDistribution<RealType>::operator()(Random& r, RealType mu,
85  RealType sigma) const throw() {
86  // N.B. These constants can be regarded as "exact", so that the same number
87  // of significant figures are used in all versions. (They serve to
88  // "bracket" the real boundary specified by the log expression.)
89  const RealType
90  m = RealType( 1.7156 ), // sqrt(8/e) (rounded up)
91  s = RealType( 0.449871), // Constants from Leva
92  t = RealType(-0.386595),
93  a = RealType( 0.19600 ),
94  b = RealType( 0.25472 ),
95  r1 = RealType( 0.27597 ),
96  r2 = RealType( 0.27846 );
97  RealType u, v, Q;
98  do { // This loop is executed 1.369 times on average
99  // Pick point P = (u, v)
100  u = r.template FixedU<RealType>(); // Sample u in (0,1]
101  v = m * r.template FixedS<RealType>(); // Sample v in (-m/2, m/2); avoid 0
102  // Compute quadratic form Q
103  const RealType x = u - s;
104  const RealType y = (v < 0 ? -v : v) - t; // Sun has no long double abs!
105  Q = x*x + y * (a*y - b*x);
106  } while ( Q >= r1 && // accept P if Q < r1
107  ( Q > r2 || // reject P if Q > r2
108  v*v > - 4 * u*u * std::log(u) ) ); // accept P if v^2 <= ...
109  return mu + sigma * (v / u); // return the slope of P (note u != 0)
110  }
111 
112 } // namespace RandomLib
113 
114 #endif // RANDOMLIB_NORMALDISTRIBUTION_HPP
Generate random integers, reals, and booleans.
RealType operator()(Random &r, RealType mu=RealType(0), RealType sigma=RealType(1)) const