1/* boost random/normal_distribution.hpp header file
2 *
3 * Copyright Jens Maurer 2000-2001
4 * Copyright Steven Watanabe 2010-2011
5 * Distributed under the Boost Software License, Version 1.0. (See
6 * accompanying file LICENSE_1_0.txt or copy at
7 * http://www.boost.org/LICENSE_1_0.txt)
8 *
9 * See http://www.boost.org for most recent version including documentation.
10 *
11 * $Id$
12 *
13 * Revision history
14 * 2001-02-18 moved to individual header files
15 */
16
17#ifndef BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP
18#define BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP
19
20#include <boost/config/no_tr1/cmath.hpp>
21#include <istream>
22#include <iosfwd>
23#include <boost/assert.hpp>
24#include <boost/limits.hpp>
25#include <boost/static_assert.hpp>
26#include <boost/random/detail/config.hpp>
27#include <boost/random/detail/operators.hpp>
28#include <boost/random/detail/int_float_pair.hpp>
29#include <boost/random/uniform_01.hpp>
30#include <boost/random/exponential_distribution.hpp>
31
32namespace boost {
33namespace random {
34
35namespace detail {
36
37// tables for the ziggurat algorithm
38template<class RealType>
39struct normal_table {
40 static const RealType table_x[129];
41 static const RealType table_y[129];
42};
43
44template<class RealType>
45const RealType normal_table<RealType>::table_x[129] = {
46 3.7130862467403632609, 3.4426198558966521214, 3.2230849845786185446, 3.0832288582142137009,
47 2.9786962526450169606, 2.8943440070186706210, 2.8231253505459664379, 2.7611693723841538514,
48 2.7061135731187223371, 2.6564064112581924999, 2.6109722484286132035, 2.5690336259216391328,
49 2.5300096723854666170, 2.4934545220919507609, 2.4590181774083500943, 2.4264206455302115930,
50 2.3954342780074673425, 2.3658713701139875435, 2.3375752413355307354, 2.3104136836950021558,
51 2.2842740596736568056, 2.2590595738653295251, 2.2346863955870569803, 2.2110814088747278106,
52 2.1881804320720206093, 2.1659267937448407377, 2.1442701823562613518, 2.1231657086697899595,
53 2.1025731351849988838, 2.0824562379877246441, 2.0627822745039633575, 2.0435215366506694976,
54 2.0246469733729338782, 2.0061338699589668403, 1.9879595741230607243, 1.9701032608497132242,
55 1.9525457295488889058, 1.9352692282919002011, 1.9182573008597320303, 1.9014946531003176140,
56 1.8849670357028692380, 1.8686611409895420085, 1.8525645117230870617, 1.8366654602533840447,
57 1.8209529965910050740, 1.8054167642140487420, 1.7900469825946189862, 1.7748343955807692457,
58 1.7597702248942318749, 1.7448461281083765085, 1.7300541605582435350, 1.7153867407081165482,
59 1.7008366185643009437, 1.6863968467734863258, 1.6720607540918522072, 1.6578219209482075462,
60 1.6436741568569826489, 1.6296114794646783962, 1.6156280950371329644, 1.6017183802152770587,
61 1.5878768648844007019, 1.5740982160167497219, 1.5603772223598406870, 1.5467087798535034608,
62 1.5330878776675560787, 1.5195095847593707806, 1.5059690368565502602, 1.4924614237746154081,
63 1.4789819769830978546, 1.4655259573357946276, 1.4520886428822164926, 1.4386653166774613138,
64 1.4252512545068615734, 1.4118417124397602509, 1.3984319141236063517, 1.3850170377251486449,
65 1.3715922024197322698, 1.3581524543224228739, 1.3446927517457130432, 1.3312079496576765017,
66 1.3176927832013429910, 1.3041418501204215390, 1.2905495919178731508, 1.2769102735516997175,
67 1.2632179614460282310, 1.2494664995643337480, 1.2356494832544811749, 1.2217602305309625678,
68 1.2077917504067576028, 1.1937367078237721994, 1.1795873846544607035, 1.1653356361550469083,
69 1.1509728421389760651, 1.1364898520030755352, 1.1218769225722540661, 1.1071236475235353980,
70 1.0922188768965537614, 1.0771506248819376573, 1.0619059636836193998, 1.0464709007525802629,
71 1.0308302360564555907, 1.0149673952392994716, 0.99886423348064351303, 0.98250080350276038481,
72 0.96585507938813059489, 0.94890262549791195381, 0.93161619660135381056, 0.91396525100880177644,
73 0.89591535256623852894, 0.87742742909771569142, 0.85845684317805086354, 0.83895221428120745572,
74 0.81885390668331772331, 0.79809206062627480454, 0.77658398787614838598, 0.75423066443451007146,
75 0.73091191062188128150, 0.70647961131360803456, 0.68074791864590421664, 0.65347863871504238702,
76 0.62435859730908822111, 0.59296294244197797913, 0.55869217837551797140, 0.52065603872514491759,
77 0.47743783725378787681, 0.42654798630330512490, 0.36287143102841830424, 0.27232086470466385065,
78 0
79};
80
81template<class RealType>
82const RealType normal_table<RealType>::table_y[129] = {
83 0, 0.0026696290839025035092, 0.0055489952208164705392, 0.0086244844129304709682,
84 0.011839478657982313715, 0.015167298010672042468, 0.018592102737165812650, 0.022103304616111592615,
85 0.025693291936149616572, 0.029356317440253829618, 0.033087886146505155566, 0.036884388786968774128,
86 0.040742868074790604632, 0.044660862200872429800, 0.048636295860284051878, 0.052667401903503169793,
87 0.056752663481538584188, 0.060890770348566375972, 0.065080585213631873753, 0.069321117394180252601,
88 0.073611501884754893389, 0.077950982514654714188, 0.082338898242957408243, 0.086774671895542968998,
89 0.091257800827634710201, 0.09578784912257815216, 0.10036444102954554013, 0.10498725541035453978,
90 0.10965602101581776100, 0.11437051244988827452, 0.11913054670871858767, 0.12393598020398174246,
91 0.12878670619710396109, 0.13368265258464764118, 0.13862377998585103702, 0.14361008009193299469,
92 0.14864157424369696566, 0.15371831220958657066, 0.15884037114093507813, 0.16400785468492774791,
93 0.16922089223892475176, 0.17447963833240232295, 0.17978427212496211424, 0.18513499701071343216,
94 0.19053204032091372112, 0.19597565311811041399, 0.20146611007620324118, 0.20700370944187380064,
95 0.21258877307373610060, 0.21822164655637059599, 0.22390269938713388747, 0.22963232523430270355,
96 0.23541094226572765600, 0.24123899354775131610, 0.24711694751469673582, 0.25304529850976585934,
97 0.25902456739871074263, 0.26505530225816194029, 0.27113807914102527343, 0.27727350292189771153,
98 0.28346220822601251779, 0.28970486044581049771, 0.29600215684985583659, 0.30235482778947976274,
99 0.30876363800925192282, 0.31522938806815752222, 0.32175291587920862031, 0.32833509837615239609,
100 0.33497685331697116147, 0.34167914123501368412, 0.34844296754987246935, 0.35526938485154714435,
101 0.36215949537303321162, 0.36911445366827513952, 0.37613546951445442947, 0.38322381105988364587,
102 0.39038080824138948916, 0.39760785649804255208, 0.40490642081148835099, 0.41227804010702462062,
103 0.41972433205403823467, 0.42724699830956239880, 0.43484783025466189638, 0.44252871528024661483,
104 0.45029164368692696086, 0.45813871627287196483, 0.46607215269457097924, 0.47409430069824960453,
105 0.48220764633483869062, 0.49041482528932163741, 0.49871863547658432422, 0.50712205108130458951,
106 0.51562823824987205196, 0.52424057267899279809, 0.53296265938998758838, 0.54179835503172412311,
107 0.55075179312105527738, 0.55982741271069481791, 0.56902999107472161225, 0.57836468112670231279,
108 0.58783705444182052571, 0.59745315095181228217, 0.60721953663260488551, 0.61714337082656248870,
109 0.62723248525781456578, 0.63749547734314487428, 0.64794182111855080873, 0.65858200005865368016,
110 0.66942766735770616891, 0.68049184100641433355, 0.69178914344603585279, 0.70333609902581741633,
111 0.71515150742047704368, 0.72725691835450587793, 0.73967724368333814856, 0.75244155918570380145,
112 0.76558417390923599480, 0.77914608594170316563, 0.79317701178385921053, 0.80773829469612111340,
113 0.82290721139526200050, 0.83878360531064722379, 0.85550060788506428418, 0.87324304892685358879,
114 0.89228165080230272301, 0.91304364799203805999, 0.93628268170837107547, 0.96359969315576759960,
115 1
116};
117
118
119template<class RealType = double>
120struct unit_normal_distribution
121{
122 template<class Engine>
123 RealType operator()(Engine& eng) {
124 const double * const table_x = normal_table<double>::table_x;
125 const double * const table_y = normal_table<double>::table_y;
126 for(;;) {
127 std::pair<RealType, int> vals = generate_int_float_pair<RealType, 8>(eng);
128 int i = vals.second;
129 int sign = (i & 1) * 2 - 1;
130 i = i >> 1;
131 RealType x = vals.first * RealType(table_x[i]);
132 if(x < table_x[i + 1]) return x * sign;
133 if(i == 0) return generate_tail(eng) * sign;
134
135 RealType y01 = uniform_01<RealType>()(eng);
136 RealType y = RealType(table_y[i]) + y01 * RealType(table_y[i + 1] - table_y[i]);
137
138 // These store the value y - bound, or something proportional to that difference:
139 RealType y_above_ubound, y_above_lbound;
140
141 // There are three cases to consider:
142 // - convex regions (where x[i] > x[j] >= 1)
143 // - concave regions (where 1 <= x[i] < x[j])
144 // - region containing the inflection point (where x[i] > 1 > x[j])
145 // For convex (concave), exp^(-x^2/2) is bounded below (above) by the tangent at
146 // (x[i],y[i]) and is bounded above (below) by the diagonal line from (x[i+1],y[i+1]) to
147 // (x[i],y[i]).
148 //
149 // *If* the inflection point region satisfies slope(x[i+1]) < slope(diagonal), then we
150 // can treat the inflection region as a convex region: this condition is necessary and
151 // sufficient to ensure that the curve lies entirely below the diagonal (that, in turn,
152 // also implies that it will be above the tangent at x[i]).
153 //
154 // For the current table size (128), this is satisfied: slope(x[i+1]) = -0.60653 <
155 // slope(diag) = -0.60649, and so we have only two cases below instead of three.
156
157 if (table_x[i] >= 1) { // convex (incl. inflection)
158 y_above_ubound = RealType(table_x[i] - table_x[i+1]) * y01 - (RealType(table_x[i]) - x);
159 y_above_lbound = y - (RealType(table_y[i]) + (RealType(table_x[i]) - x) * RealType(table_y[i]) * RealType(table_x[i]));
160 }
161 else { // concave
162 y_above_lbound = RealType(table_x[i] - table_x[i+1]) * y01 - (RealType(table_x[i]) - x);
163 y_above_ubound = y - (RealType(table_y[i]) + (RealType(table_x[i]) - x) * RealType(table_y[i]) * RealType(table_x[i]));
164 }
165
166 if (y_above_ubound < 0 // if above the upper bound reject immediately
167 &&
168 (
169 y_above_lbound < 0 // If below the lower bound accept immediately
170 ||
171 y < f(x) // Otherwise it's between the bounds and we need a full check
172 )
173 ) {
174 return x * sign;
175 }
176 }
177 }
178 static RealType f(RealType x) {
179 using std::exp;
180 return exp(-(x*x/2));
181 }
182 // Generate from the tail using rejection sampling from the exponential(x_1) distribution,
183 // shifted by x_1. This looks a little different from the usual rejection sampling because it
184 // transforms the condition by taking the log of both sides, thus avoiding the costly exp() call
185 // on the RHS, then takes advantage of the fact that -log(unif01) is simply generating an
186 // exponential (by inverse cdf sampling) by replacing the log(unif01) on the LHS with a
187 // exponential(1) draw, y.
188 template<class Engine>
189 RealType generate_tail(Engine& eng) {
190 const RealType tail_start = RealType(normal_table<double>::table_x[1]);
191 boost::random::exponential_distribution<RealType> exp_x(tail_start);
192 boost::random::exponential_distribution<RealType> exp_y;
193 for(;;) {
194 RealType x = exp_x(eng);
195 RealType y = exp_y(eng);
196 // If we were doing non-transformed rejection sampling, this condition would be:
197 // if (unif01 < exp(-.5*x*x)) return x + tail_start;
198 if(2*y > x*x) return x + tail_start;
199 }
200 }
201};
202
203} // namespace detail
204
205
206/**
207 * Instantiations of class template normal_distribution model a
208 * \random_distribution. Such a distribution produces random numbers
209 * @c x distributed with probability density function
210 * \f$\displaystyle p(x) =
211 * \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(x-\mu)^2}{2\sigma^2}}
212 * \f$,
213 * where mean and sigma are the parameters of the distribution.
214 *
215 * The implementation uses the "ziggurat" algorithm, as described in
216 *
217 * @blockquote
218 * "The Ziggurat Method for Generating Random Variables",
219 * George Marsaglia and Wai Wan Tsang, Journal of Statistical Software,
220 * Volume 5, Number 8 (2000), 1-7.
221 * @endblockquote
222 */
223template<class RealType = double>
224class normal_distribution
225{
226public:
227 typedef RealType input_type;
228 typedef RealType result_type;
229
230 class param_type {
231 public:
232 typedef normal_distribution distribution_type;
233
234 /**
235 * Constructs a @c param_type with a given mean and
236 * standard deviation.
237 *
238 * Requires: sigma >= 0
239 */
240 explicit param_type(RealType mean_arg = RealType(0.0),
241 RealType sigma_arg = RealType(1.0))
242 : _mean(mean_arg),
243 _sigma(sigma_arg)
244 {}
245
246 /** Returns the mean of the distribution. */
247 RealType mean() const { return _mean; }
248
249 /** Returns the standand deviation of the distribution. */
250 RealType sigma() const { return _sigma; }
251
252 /** Writes a @c param_type to a @c std::ostream. */
253 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
254 { os << parm._mean << " " << parm._sigma ; return os; }
255
256 /** Reads a @c param_type from a @c std::istream. */
257 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
258 { is >> parm._mean >> std::ws >> parm._sigma; return is; }
259
260 /** Returns true if the two sets of parameters are the same. */
261 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
262 { return lhs._mean == rhs._mean && lhs._sigma == rhs._sigma; }
263
264 /** Returns true if the two sets of parameters are the different. */
265 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
266
267 private:
268 RealType _mean;
269 RealType _sigma;
270 };
271
272 /**
273 * Constructs a @c normal_distribution object. @c mean and @c sigma are
274 * the parameters for the distribution.
275 *
276 * Requires: sigma >= 0
277 */
278 explicit normal_distribution(const RealType& mean_arg = RealType(0.0),
279 const RealType& sigma_arg = RealType(1.0))
280 : _mean(mean_arg), _sigma(sigma_arg)
281 {
282 BOOST_ASSERT(_sigma >= RealType(0));
283 }
284
285 /**
286 * Constructs a @c normal_distribution object from its parameters.
287 */
288 explicit normal_distribution(const param_type& parm)
289 : _mean(parm.mean()), _sigma(parm.sigma())
290 {}
291
292 /** Returns the mean of the distribution. */
293 RealType mean() const { return _mean; }
294 /** Returns the standard deviation of the distribution. */
295 RealType sigma() const { return _sigma; }
296
297 /** Returns the smallest value that the distribution can produce. */
298 RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const
299 { return -std::numeric_limits<RealType>::infinity(); }
300 /** Returns the largest value that the distribution can produce. */
301 RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
302 { return std::numeric_limits<RealType>::infinity(); }
303
304 /** Returns the parameters of the distribution. */
305 param_type param() const { return param_type(_mean, _sigma); }
306 /** Sets the parameters of the distribution. */
307 void param(const param_type& parm)
308 {
309 _mean = parm.mean();
310 _sigma = parm.sigma();
311 }
312
313 /**
314 * Effects: Subsequent uses of the distribution do not depend
315 * on values produced by any engine prior to invoking reset.
316 */
317 void reset() { }
318
319 /** Returns a normal variate. */
320 template<class Engine>
321 result_type operator()(Engine& eng)
322 {
323 detail::unit_normal_distribution<RealType> impl;
324 return impl(eng) * _sigma + _mean;
325 }
326
327 /** Returns a normal variate with parameters specified by @c param. */
328 template<class URNG>
329 result_type operator()(URNG& urng, const param_type& parm)
330 {
331 return normal_distribution(parm)(urng);
332 }
333
334 /** Writes a @c normal_distribution to a @c std::ostream. */
335 BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, normal_distribution, nd)
336 {
337 os << nd._mean << " " << nd._sigma;
338 return os;
339 }
340
341 /** Reads a @c normal_distribution from a @c std::istream. */
342 BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, normal_distribution, nd)
343 {
344 is >> std::ws >> nd._mean >> std::ws >> nd._sigma;
345 return is;
346 }
347
348 /**
349 * Returns true if the two instances of @c normal_distribution will
350 * return identical sequences of values given equal generators.
351 */
352 BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(normal_distribution, lhs, rhs)
353 {
354 return lhs._mean == rhs._mean && lhs._sigma == rhs._sigma;
355 }
356
357 /**
358 * Returns true if the two instances of @c normal_distribution will
359 * return different sequences of values given equal generators.
360 */
361 BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(normal_distribution)
362
363private:
364 RealType _mean, _sigma;
365
366};
367
368} // namespace random
369
370using random::normal_distribution;
371
372} // namespace boost
373
374#endif // BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP
375