boost/random/normal_distribution.hpp

0001 /* boost random/normal_distribution.hpp header file
0002  *
0003  * Copyright Jens Maurer 2000-2001
0004  * Copyright Steven Watanabe 2010-2011
0005  * Distributed under the Boost Software License, Version 1.0. (See
0006  * accompanying file LICENSE_1_0.txt or copy at
0007  * http://www.boost.org/LICENSE_1_0.txt)
0008  *
0009  * See http://www.boost.org for most recent version including documentation.
0010  *
0011  * $Id$
0012  *
0013  * Revision history
0014  *  2001-02-18  moved to individual header files
0015  */
0016
0017 #ifndef BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP
0018 #define BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP
0019
0020 #include <boost/config/no_tr1/cmath.hpp>
0021 #include <istream>
0022 #include <iosfwd>
0023 #include <boost/assert.hpp>
0024 #include <boost/limits.hpp>
0025 #include <boost/static_assert.hpp>
0026 #include <boost/random/detail/config.hpp>
0027 #include <boost/random/detail/operators.hpp>
0028 #include <boost/random/detail/int_float_pair.hpp>
0029 #include <boost/random/uniform_01.hpp>
0030 #include <boost/random/exponential_distribution.hpp>
0031
0032 namespace boost {
0033 namespace random {
0034
0035 namespace detail {
0036
0037 // tables for the ziggurat algorithm
0038 template<class RealType>
0039 struct normal_table {
0040     static const RealType table_x[129];
0041     static const RealType table_y[129];
0042 };
0043
0044 template<class RealType>
0045 const RealType normal_table<RealType>::table_x[129] = {
0046     3.7130862467403632609, 3.4426198558966521214, 3.2230849845786185446, 3.0832288582142137009,
0047     2.9786962526450169606, 2.8943440070186706210, 2.8231253505459664379, 2.7611693723841538514,
0048     2.7061135731187223371, 2.6564064112581924999, 2.6109722484286132035, 2.5690336259216391328,
0049     2.5300096723854666170, 2.4934545220919507609, 2.4590181774083500943, 2.4264206455302115930,
0050     2.3954342780074673425, 2.3658713701139875435, 2.3375752413355307354, 2.3104136836950021558,
0051     2.2842740596736568056, 2.2590595738653295251, 2.2346863955870569803, 2.2110814088747278106,
0052     2.1881804320720206093, 2.1659267937448407377, 2.1442701823562613518, 2.1231657086697899595,
0053     2.1025731351849988838, 2.0824562379877246441, 2.0627822745039633575, 2.0435215366506694976,
0054     2.0246469733729338782, 2.0061338699589668403, 1.9879595741230607243, 1.9701032608497132242,
0055     1.9525457295488889058, 1.9352692282919002011, 1.9182573008597320303, 1.9014946531003176140,
0056     1.8849670357028692380, 1.8686611409895420085, 1.8525645117230870617, 1.8366654602533840447,
0057     1.8209529965910050740, 1.8054167642140487420, 1.7900469825946189862, 1.7748343955807692457,
0058     1.7597702248942318749, 1.7448461281083765085, 1.7300541605582435350, 1.7153867407081165482,
0059     1.7008366185643009437, 1.6863968467734863258, 1.6720607540918522072, 1.6578219209482075462,
0060     1.6436741568569826489, 1.6296114794646783962, 1.6156280950371329644, 1.6017183802152770587,
0061     1.5878768648844007019, 1.5740982160167497219, 1.5603772223598406870, 1.5467087798535034608,
0062     1.5330878776675560787, 1.5195095847593707806, 1.5059690368565502602, 1.4924614237746154081,
0063     1.4789819769830978546, 1.4655259573357946276, 1.4520886428822164926, 1.4386653166774613138,
0064     1.4252512545068615734, 1.4118417124397602509, 1.3984319141236063517, 1.3850170377251486449,
0065     1.3715922024197322698, 1.3581524543224228739, 1.3446927517457130432, 1.3312079496576765017,
0066     1.3176927832013429910, 1.3041418501204215390, 1.2905495919178731508, 1.2769102735516997175,
0067     1.2632179614460282310, 1.2494664995643337480, 1.2356494832544811749, 1.2217602305309625678,
0068     1.2077917504067576028, 1.1937367078237721994, 1.1795873846544607035, 1.1653356361550469083,
0069     1.1509728421389760651, 1.1364898520030755352, 1.1218769225722540661, 1.1071236475235353980,
0070     1.0922188768965537614, 1.0771506248819376573, 1.0619059636836193998, 1.0464709007525802629,
0071     1.0308302360564555907, 1.0149673952392994716, 0.99886423348064351303, 0.98250080350276038481,
0072     0.96585507938813059489, 0.94890262549791195381, 0.93161619660135381056, 0.91396525100880177644,
0073     0.89591535256623852894, 0.87742742909771569142, 0.85845684317805086354, 0.83895221428120745572,
0074     0.81885390668331772331, 0.79809206062627480454, 0.77658398787614838598, 0.75423066443451007146,
0075     0.73091191062188128150, 0.70647961131360803456, 0.68074791864590421664, 0.65347863871504238702,
0076     0.62435859730908822111, 0.59296294244197797913, 0.55869217837551797140, 0.52065603872514491759,
0077     0.47743783725378787681, 0.42654798630330512490, 0.36287143102841830424, 0.27232086470466385065,
0078     0
0079 };
0080
0081 template<class RealType>
0082 const RealType normal_table<RealType>::table_y[129] = {
0083     0, 0.0026696290839025035092, 0.0055489952208164705392, 0.0086244844129304709682,
0084     0.011839478657982313715, 0.015167298010672042468, 0.018592102737165812650, 0.022103304616111592615,
0085     0.025693291936149616572, 0.029356317440253829618, 0.033087886146505155566, 0.036884388786968774128,
0086     0.040742868074790604632, 0.044660862200872429800, 0.048636295860284051878, 0.052667401903503169793,
0087     0.056752663481538584188, 0.060890770348566375972, 0.065080585213631873753, 0.069321117394180252601,
0088     0.073611501884754893389, 0.077950982514654714188, 0.082338898242957408243, 0.086774671895542968998,
0089     0.091257800827634710201, 0.09578784912257815216, 0.10036444102954554013, 0.10498725541035453978,
0090     0.10965602101581776100, 0.11437051244988827452, 0.11913054670871858767, 0.12393598020398174246,
0091     0.12878670619710396109, 0.13368265258464764118, 0.13862377998585103702, 0.14361008009193299469,
0092     0.14864157424369696566, 0.15371831220958657066, 0.15884037114093507813, 0.16400785468492774791,
0093     0.16922089223892475176, 0.17447963833240232295, 0.17978427212496211424, 0.18513499701071343216,
0094     0.19053204032091372112, 0.19597565311811041399, 0.20146611007620324118, 0.20700370944187380064,
0095     0.21258877307373610060, 0.21822164655637059599, 0.22390269938713388747, 0.22963232523430270355,
0096     0.23541094226572765600, 0.24123899354775131610, 0.24711694751469673582, 0.25304529850976585934,
0097     0.25902456739871074263, 0.26505530225816194029, 0.27113807914102527343, 0.27727350292189771153,
0098     0.28346220822601251779, 0.28970486044581049771, 0.29600215684985583659, 0.30235482778947976274,
0099     0.30876363800925192282, 0.31522938806815752222, 0.32175291587920862031, 0.32833509837615239609,
0100     0.33497685331697116147, 0.34167914123501368412, 0.34844296754987246935, 0.35526938485154714435,
0101     0.36215949537303321162, 0.36911445366827513952, 0.37613546951445442947, 0.38322381105988364587,
0102     0.39038080824138948916, 0.39760785649804255208, 0.40490642081148835099, 0.41227804010702462062,
0103     0.41972433205403823467, 0.42724699830956239880, 0.43484783025466189638, 0.44252871528024661483,
0104     0.45029164368692696086, 0.45813871627287196483, 0.46607215269457097924, 0.47409430069824960453,
0105     0.48220764633483869062, 0.49041482528932163741, 0.49871863547658432422, 0.50712205108130458951,
0106     0.51562823824987205196, 0.52424057267899279809, 0.53296265938998758838, 0.54179835503172412311,
0107     0.55075179312105527738, 0.55982741271069481791, 0.56902999107472161225, 0.57836468112670231279,
0108     0.58783705444182052571, 0.59745315095181228217, 0.60721953663260488551, 0.61714337082656248870,
0109     0.62723248525781456578, 0.63749547734314487428, 0.64794182111855080873, 0.65858200005865368016,
0110     0.66942766735770616891, 0.68049184100641433355, 0.69178914344603585279, 0.70333609902581741633,
0111     0.71515150742047704368, 0.72725691835450587793, 0.73967724368333814856, 0.75244155918570380145,
0112     0.76558417390923599480, 0.77914608594170316563, 0.79317701178385921053, 0.80773829469612111340,
0113     0.82290721139526200050, 0.83878360531064722379, 0.85550060788506428418, 0.87324304892685358879,
0114     0.89228165080230272301, 0.91304364799203805999, 0.93628268170837107547, 0.96359969315576759960,
0115     1
0116 };
0117
0118
0119 template<class RealType = double>
0120 struct unit_normal_distribution
0121 {
0122     template<class Engine>
0123     RealType operator()(Engine& eng) {
0124         const double * const table_x = normal_table<double>::table_x;
0125         const double * const table_y = normal_table<double>::table_y;
0126         for(;;) {
0127             std::pair<RealType, int> vals = generate_int_float_pair<RealType, 8>(eng);
0128             int i = vals.second;
0129             int sign = (i & 1) * 2 - 1;
0130             i = i >> 1;
0131             RealType x = vals.first * RealType(table_x[i]);
0132             if(x < table_x[i + 1]) return x * sign;
0133             if(i == 0) return generate_tail(eng) * sign;
0134
0135             RealType y01 = uniform_01<RealType>()(eng);
0136             RealType y = RealType(table_y[i]) + y01 * RealType(table_y[i + 1] - table_y[i]);
0137
0138             // These store the value y - bound, or something proportional to that difference:
0139             RealType y_above_ubound, y_above_lbound;
0140
0141             // There are three cases to consider:
0142             // - convex regions (where x[i] > x[j] >= 1)
0143             // - concave regions (where 1 <= x[i] < x[j])
0144             // - region containing the inflection point (where x[i] > 1 > x[j])
0145             // For convex (concave), exp^(-x^2/2) is bounded below (above) by the tangent at
0146             // (x[i],y[i]) and is bounded above (below) by the diagonal line from (x[i+1],y[i+1]) to
0147             // (x[i],y[i]).
0148             //
0149             // *If* the inflection point region satisfies slope(x[i+1]) < slope(diagonal), then we
0150             // can treat the inflection region as a convex region: this condition is necessary and
0151             // sufficient to ensure that the curve lies entirely below the diagonal (that, in turn,
0152             // also implies that it will be above the tangent at x[i]).
0153             //
0154             // For the current table size (128), this is satisfied: slope(x[i+1]) = -0.60653 <
0155             // slope(diag) = -0.60649, and so we have only two cases below instead of three.
0156
0157             if (table_x[i] >= 1) { // convex (incl. inflection)
0158                 y_above_ubound = RealType(table_x[i] - table_x[i+1]) * y01 - (RealType(table_x[i]) - x);
0159                 y_above_lbound = y - (RealType(table_y[i]) + (RealType(table_x[i]) - x) * RealType(table_y[i]) * RealType(table_x[i]));
0160             }
0161             else { // concave
0162                 y_above_lbound = RealType(table_x[i] - table_x[i+1]) * y01 - (RealType(table_x[i]) - x);
0163                 y_above_ubound = y - (RealType(table_y[i]) + (RealType(table_x[i]) - x) * RealType(table_y[i]) * RealType(table_x[i]));
0164             }
0165
0166             if (y_above_ubound < 0 // if above the upper bound reject immediately
0167                     &&
0168                     (
0169                      y_above_lbound < 0 // If below the lower bound accept immediately
0170                      ||
0171                      y < f(x) // Otherwise it's between the bounds and we need a full check
0172                     )
0173                ) {
0174                 return x * sign;
0175             }
0176         }
0177     }
0178     static RealType f(RealType x) {
0179         using std::exp;
0180         return exp(-(x*x/2));
0181     }
0182     // Generate from the tail using rejection sampling from the exponential(x_1) distribution,
0183     // shifted by x_1.  This looks a little different from the usual rejection sampling because it
0184     // transforms the condition by taking the log of both sides, thus avoiding the costly exp() call
0185     // on the RHS, then takes advantage of the fact that -log(unif01) is simply generating an
0186     // exponential (by inverse cdf sampling) by replacing the log(unif01) on the LHS with a
0187     // exponential(1) draw, y.
0188     template<class Engine>
0189     RealType generate_tail(Engine& eng) {
0190         const RealType tail_start = RealType(normal_table<double>::table_x[1]);
0191         boost::random::exponential_distribution<RealType> exp_x(tail_start);
0192         boost::random::exponential_distribution<RealType> exp_y;
0193         for(;;) {
0194             RealType x = exp_x(eng);
0195             RealType y = exp_y(eng);
0196             // If we were doing non-transformed rejection sampling, this condition would be:
0197             // if (unif01 < exp(-.5*x*x)) return x + tail_start;
0198             if(2*y > x*x) return x + tail_start;
0199         }
0200     }
0201 };
0202
0203 } // namespace detail
0204
0205
0206 /**
0207  * Instantiations of class template normal_distribution model a
0208  * \random_distribution. Such a distribution produces random numbers
0209  * @c x distributed with probability density function
0210  * \f$\displaystyle p(x) =
0211  *   \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(x-\mu)^2}{2\sigma^2}}
0212  * \f$,
0213  * where mean and sigma are the parameters of the distribution.
0214  *
0215  * The implementation uses the "ziggurat" algorithm, as described in
0216  *
0217  *  @blockquote
0218  *  "The Ziggurat Method for Generating Random Variables",
0219  *  George Marsaglia and Wai Wan Tsang, Journal of Statistical Software,
0220  *  Volume 5, Number 8 (2000), 1-7.
0221  *  @endblockquote
0222  */
0223 template<class RealType = double>
0224 class normal_distribution
0225 {
0226 public:
0227     typedef RealType input_type;
0228     typedef RealType result_type;
0229
0230     class param_type {
0231     public:
0232         typedef normal_distribution distribution_type;
0233
0234         /**
0235          * Constructs a @c param_type with a given mean and
0236          * standard deviation.
0237          *
0238          * Requires: sigma >= 0
0239          */
0240         explicit param_type(RealType mean_arg = RealType(0.0),
0241                             RealType sigma_arg = RealType(1.0))
0242           : _mean(mean_arg),
0243             _sigma(sigma_arg)
0244         {}
0245
0246         /** Returns the mean of the distribution. */
0247         RealType mean() const { return _mean; }
0248
0249         /** Returns the standand deviation of the distribution. */
0250         RealType sigma() const { return _sigma; }
0251
0252         /** Writes a @c param_type to a @c std::ostream. */
0253         BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
0254         { os << parm._mean << " " << parm._sigma ; return os; }
0255
0256         /** Reads a @c param_type from a @c std::istream. */
0257         BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
0258         { is >> parm._mean >> std::ws >> parm._sigma; return is; }
0259
0260         /** Returns true if the two sets of parameters are the same. */
0261         BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
0262         { return lhs._mean == rhs._mean && lhs._sigma == rhs._sigma; }
0263
0264         /** Returns true if the two sets of parameters are the different. */
0265         BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
0266
0267     private:
0268         RealType _mean;
0269         RealType _sigma;
0270     };
0271
0272     /**
0273      * Constructs a @c normal_distribution object. @c mean and @c sigma are
0274      * the parameters for the distribution.
0275      *
0276      * Requires: sigma >= 0
0277      */
0278     explicit normal_distribution(const RealType& mean_arg = RealType(0.0),
0279                                  const RealType& sigma_arg = RealType(1.0))
0280       : _mean(mean_arg), _sigma(sigma_arg)
0281     {
0282         BOOST_ASSERT(_sigma >= RealType(0));
0283     }
0284
0285     /**
0286      * Constructs a @c normal_distribution object from its parameters.
0287      */
0288     explicit normal_distribution(const param_type& parm)
0289       : _mean(parm.mean()), _sigma(parm.sigma())
0290     {}
0291
0292     /**  Returns the mean of the distribution. */
0293     RealType mean() const { return _mean; }
0294     /** Returns the standard deviation of the distribution. */
0295     RealType sigma() const { return _sigma; }
0296
0297     /** Returns the smallest value that the distribution can produce. */
0298     RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const
0299     { return -std::numeric_limits<RealType>::infinity(); }
0300     /** Returns the largest value that the distribution can produce. */
0301     RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
0302     { return std::numeric_limits<RealType>::infinity(); }
0303
0304     /** Returns the parameters of the distribution. */
0305     param_type param() const { return param_type(_mean, _sigma); }
0306     /** Sets the parameters of the distribution. */
0307     void param(const param_type& parm)
0308     {
0309         _mean = parm.mean();
0310         _sigma = parm.sigma();
0311     }
0312
0313     /**
0314      * Effects: Subsequent uses of the distribution do not depend
0315      * on values produced by any engine prior to invoking reset.
0316      */
0317     void reset() { }
0318
0319     /**  Returns a normal variate. */
0320     template<class Engine>
0321     result_type operator()(Engine& eng)
0322     {
0323         detail::unit_normal_distribution<RealType> impl;
0324         return impl(eng) * _sigma + _mean;
0325     }
0326
0327     /** Returns a normal variate with parameters specified by @c param. */
0328     template<class URNG>
0329     result_type operator()(URNG& urng, const param_type& parm)
0330     {
0331         return normal_distribution(parm)(urng);
0332     }
0333
0334     /** Writes a @c normal_distribution to a @c std::ostream. */
0335     BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, normal_distribution, nd)
0336     {
0337         os << nd._mean << " " << nd._sigma;
0338         return os;
0339     }
0340
0341     /** Reads a @c normal_distribution from a @c std::istream. */
0342     BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, normal_distribution, nd)
0343     {
0344         is >> std::ws >> nd._mean >> std::ws >> nd._sigma;
0345         return is;
0346     }
0347
0348     /**
0349      * Returns true if the two instances of @c normal_distribution will
0350      * return identical sequences of values given equal generators.
0351      */
0352     BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(normal_distribution, lhs, rhs)
0353     {
0354         return lhs._mean == rhs._mean && lhs._sigma == rhs._sigma;
0355     }
0356
0357     /**
0358      * Returns true if the two instances of @c normal_distribution will
0359      * return different sequences of values given equal generators.
0360      */
0361     BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(normal_distribution)
0362
0363 private:
0364     RealType _mean, _sigma;
0365
0366 };
0367
0368 } // namespace random
0369
0370 using random::normal_distribution;
0371
0372 } // namespace boost
0373
0374 #endif // BOOST_RANDOM_NORMAL_DISTRIBUTION_HPP