boost/random/gamma_distribution.hpp

   1 /* boost random/gamma_distribution.hpp header file
   2  *
   3  * Copyright Jens Maurer 2002
   4  * Copyright Steven Watanabe 2010
   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: gamma_distribution.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $
  12  *
  13  */
  14
  15 #ifndef BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
  16 #define BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
  17
  18 #include <boost/config/no_tr1/cmath.hpp>
  19 #include <istream>
  20 #include <iosfwd>
  21 #include <boost/assert.hpp>
  22 #include <boost/limits.hpp>
  23 #include <boost/static_assert.hpp>
  24 #include <boost/random/detail/config.hpp>
  25 #include <boost/random/exponential_distribution.hpp>
  26
  27 namespace boost {
  28 namespace random {
  29
  30 // The algorithm is taken from Knuth
  31
  32 /**
  33  * The gamma distribution is a continuous distribution with two
  34  * parameters alpha and beta.  It produces values > 0.
  35  *
  36  * It has
  37  * \f$\displaystyle p(x) = x^{\alpha-1}\frac{e^{-x/\beta}}{\beta^\alpha\Gamma(\alpha)}\f$.
  38  */
  39 template<class RealType = double>
  40 class gamma_distribution
  41 {
  42 public:
  43     typedef RealType input_type;
  44     typedef RealType result_type;
  45
  46     class param_type
  47     {
  48     public:
  49         typedef gamma_distribution distribution_type;
  50
  51         /**
  52          * Constructs a @c param_type object from the "alpha" and "beta"
  53          * parameters.
  54          *
  55          * Requires: alpha > 0 && beta > 0
  56          */
  57         param_type(const RealType& alpha_arg = RealType(1.0),
  58                    const RealType& beta_arg = RealType(1.0))
  59           : _alpha(alpha_arg), _beta(beta_arg)
  60         {
  61         }
  62
  63         /** Returns the "alpha" parameter of the distribution. */
  64         RealType alpha() const { return _alpha; }
  65         /** Returns the "beta" parameter of the distribution. */
  66         RealType beta() const { return _beta; }
  67
  68 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
  69         /** Writes the parameters to a @c std::ostream. */
  70         template<class CharT, class Traits>
  71         friend std::basic_ostream<CharT, Traits>&
  72         operator<<(std::basic_ostream<CharT, Traits>& os,
  73                    const param_type& parm)
  74         {
  75             os << parm._alpha << ' ' << parm._beta;
  76             return os;
  77         }
  78
  79         /** Reads the parameters from a @c std::istream. */
  80         template<class CharT, class Traits>
  81         friend std::basic_istream<CharT, Traits>&
  82         operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm)
  83         {
  84             is >> parm._alpha >> std::ws >> parm._beta;
  85             return is;
  86         }
  87 #endif
  88
  89         /** Returns true if the two sets of parameters are the same. */
  90         friend bool operator==(const param_type& lhs, const param_type& rhs)
  91         {
  92             return lhs._alpha == rhs._alpha && lhs._beta == rhs._beta;
  93         }
  94         /** Returns true if the two sets fo parameters are different. */
  95         friend bool operator!=(const param_type& lhs, const param_type& rhs)
  96         {
  97             return !(lhs == rhs);
  98         }
  99     private:
 100         RealType _alpha;
 101         RealType _beta;
 102     };
 103
 104 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
 105     BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
 106 #endif
 107
 108     /**
 109      * Creates a new gamma_distribution with parameters "alpha" and "beta".
 110      *
 111      * Requires: alpha > 0 && beta > 0
 112      */
 113     explicit gamma_distribution(const result_type& alpha_arg = result_type(1.0),
 114                                 const result_type& beta_arg = result_type(1.0))
 115       : _exp(result_type(1)), _alpha(alpha_arg), _beta(beta_arg)
 116     {
 117         BOOST_ASSERT(_alpha > result_type(0));
 118         BOOST_ASSERT(_beta > result_type(0));
 119         init();
 120     }
 121
 122     /** Constructs a @c gamma_distribution from its parameters. */
 123     explicit gamma_distribution(const param_type& parm)
 124       : _exp(result_type(1)), _alpha(parm.alpha()), _beta(parm.beta())
 125     {
 126         init();
 127     }
 128
 129     // compiler-generated copy ctor and assignment operator are fine
 130
 131     /** Returns the "alpha" paramter of the distribution. */
 132     RealType alpha() const { return _alpha; }
 133     /** Returns the "beta" parameter of the distribution. */
 134     RealType beta() const { return _beta; }
 135     /** Returns the smallest value that the distribution can produce. */
 136     RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
 137     /* Returns the largest value that the distribution can produce. */
 138     RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
 139     { return (std::numeric_limits<RealType>::infinity)(); }
 140
 141     /** Returns the parameters of the distribution. */
 142     param_type param() const { return param_type(_alpha, _beta); }
 143     /** Sets the parameters of the distribution. */
 144     void param(const param_type& parm)
 145     {
 146         _alpha = parm.alpha();
 147         _beta = parm.beta();
 148         init();
 149     }
 150
 151     /**
 152      * Effects: Subsequent uses of the distribution do not depend
 153      * on values produced by any engine prior to invoking reset.
 154      */
 155     void reset() { _exp.reset(); }
 156
 157     /**
 158      * Returns a random variate distributed according to
 159      * the gamma distribution.
 160      */
 161     template<class Engine>
 162     result_type operator()(Engine& eng)
 163     {
 164 #ifndef BOOST_NO_STDC_NAMESPACE
 165         // allow for Koenig lookup
 166         using std::tan; using std::sqrt; using std::exp; using std::log;
 167         using std::pow;
 168 #endif
 169         if(_alpha == result_type(1)) {
 170             return _exp(eng) * _beta;
 171         } else if(_alpha > result_type(1)) {
 172             // Can we have a boost::mathconst please?
 173             const result_type pi = result_type(3.14159265358979323846);
 174             for(;;) {
 175                 result_type y = tan(pi * uniform_01<RealType>()(eng));
 176                 result_type x = sqrt(result_type(2)*_alpha-result_type(1))*y
 177                     + _alpha-result_type(1);
 178                 if(x <= result_type(0))
 179                     continue;
 180                 if(uniform_01<RealType>()(eng) >
 181                     (result_type(1)+y*y) * exp((_alpha-result_type(1))
 182                                                *log(x/(_alpha-result_type(1)))
 183                                                - sqrt(result_type(2)*_alpha
 184                                                       -result_type(1))*y))
 185                     continue;
 186                 return x * _beta;
 187             }
 188         } else /* alpha < 1.0 */ {
 189             for(;;) {
 190                 result_type u = uniform_01<RealType>()(eng);
 191                 result_type y = _exp(eng);
 192                 result_type x, q;
 193                 if(u < _p) {
 194                     x = exp(-y/_alpha);
 195                     q = _p*exp(-x);
 196                 } else {
 197                     x = result_type(1)+y;
 198                     q = _p + (result_type(1)-_p) * pow(x,_alpha-result_type(1));
 199                 }
 200                 if(u >= q)
 201                     continue;
 202                 return x * _beta;
 203             }
 204         }
 205     }
 206
 207     template<class URNG>
 208     RealType operator()(URNG& urng, const param_type& parm) const
 209     {
 210         return gamma_distribution(parm)(urng);
 211     }
 212
 213 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
 214     /** Writes a @c gamma_distribution to a @c std::ostream. */
 215     template<class CharT, class Traits>
 216     friend std::basic_ostream<CharT,Traits>&
 217     operator<<(std::basic_ostream<CharT,Traits>& os,
 218                const gamma_distribution& gd)
 219     {
 220         os << gd.param();
 221         return os;
 222     }
 223
 224     /** Reads a @c gamma_distribution from a @c std::istream. */
 225     template<class CharT, class Traits>
 226     friend std::basic_istream<CharT,Traits>&
 227     operator>>(std::basic_istream<CharT,Traits>& is, gamma_distribution& gd)
 228     {
 229         gd.read(is);
 230         return is;
 231     }
 232 #endif
 233
 234     /**
 235      * Returns true if the two distributions will produce identical
 236      * sequences of random variates given equal generators.
 237      */
 238     friend bool operator==(const gamma_distribution& lhs,
 239                            const gamma_distribution& rhs)
 240     {
 241         return lhs._alpha == rhs._alpha
 242             && lhs._beta == rhs._beta
 243             && lhs._exp == rhs._exp;
 244     }
 245
 246     /**
 247      * Returns true if the two distributions can produce different
 248      * sequences of random variates, given equal generators.
 249      */
 250     friend bool operator!=(const gamma_distribution& lhs,
 251                            const gamma_distribution& rhs)
 252     {
 253         return !(lhs == rhs);
 254     }
 255
 256 private:
 257     /// \cond hide_private_members
 258
 259     template<class CharT, class Traits>
 260     void read(std::basic_istream<CharT, Traits>& is)
 261     {
 262         param_type parm;
 263         if(is >> parm) {
 264             param(parm);
 265         }
 266     }
 267
 268     void init()
 269     {
 270 #ifndef BOOST_NO_STDC_NAMESPACE
 271         // allow for Koenig lookup
 272         using std::exp;
 273 #endif
 274         _p = exp(result_type(1)) / (_alpha + exp(result_type(1)));
 275     }
 276     /// \endcond
 277
 278     exponential_distribution<RealType> _exp;
 279     result_type _alpha;
 280     result_type _beta;
 281     // some data precomputed from the parameters
 282     result_type _p;
 283 };
 284
 285
 286 } // namespace random
 287
 288 using random::gamma_distribution;
 289
 290 } // namespace boost
 291
 292 #endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP