boost/random/gamma_distribution.hpp

0001 /* boost random/gamma_distribution.hpp header file
0002  *
0003  * Copyright Jens Maurer 2002
0004  * Copyright Steven Watanabe 2010
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  */
0014
0015 #ifndef BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
0016 #define BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
0017
0018 #include <boost/config/no_tr1/cmath.hpp>
0019 #include <istream>
0020 #include <iosfwd>
0021 #include <boost/assert.hpp>
0022 #include <boost/limits.hpp>
0023 #include <boost/static_assert.hpp>
0024 #include <boost/random/detail/config.hpp>
0025 #include <boost/random/exponential_distribution.hpp>
0026
0027 namespace boost {
0028 namespace random {
0029
0030 // The algorithm is taken from Knuth
0031
0032 /**
0033  * The gamma distribution is a continuous distribution with two
0034  * parameters alpha and beta.  It produces values > 0.
0035  *
0036  * It has
0037  * \f$\displaystyle p(x) = x^{\alpha-1}\frac{e^{-x/\beta}}{\beta^\alpha\Gamma(\alpha)}\f$.
0038  */
0039 template<class RealType = double>
0040 class gamma_distribution
0041 {
0042 public:
0043     typedef RealType input_type;
0044     typedef RealType result_type;
0045
0046     class param_type
0047     {
0048     public:
0049         typedef gamma_distribution distribution_type;
0050
0051         /**
0052          * Constructs a @c param_type object from the "alpha" and "beta"
0053          * parameters.
0054          *
0055          * Requires: alpha > 0 && beta > 0
0056          */
0057         param_type(const RealType& alpha_arg = RealType(1.0),
0058                    const RealType& beta_arg = RealType(1.0))
0059           : _alpha(alpha_arg), _beta(beta_arg)
0060         {
0061         }
0062
0063         /** Returns the "alpha" parameter of the distribution. */
0064         RealType alpha() const { return _alpha; }
0065         /** Returns the "beta" parameter of the distribution. */
0066         RealType beta() const { return _beta; }
0067
0068 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
0069         /** Writes the parameters to a @c std::ostream. */
0070         template<class CharT, class Traits>
0071         friend std::basic_ostream<CharT, Traits>&
0072         operator<<(std::basic_ostream<CharT, Traits>& os,
0073                    const param_type& parm)
0074         {
0075             os << parm._alpha << ' ' << parm._beta;
0076             return os;
0077         }
0078
0079         /** Reads the parameters from a @c std::istream. */
0080         template<class CharT, class Traits>
0081         friend std::basic_istream<CharT, Traits>&
0082         operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm)
0083         {
0084             is >> parm._alpha >> std::ws >> parm._beta;
0085             return is;
0086         }
0087 #endif
0088
0089         /** Returns true if the two sets of parameters are the same. */
0090         friend bool operator==(const param_type& lhs, const param_type& rhs)
0091         {
0092             return lhs._alpha == rhs._alpha && lhs._beta == rhs._beta;
0093         }
0094         /** Returns true if the two sets fo parameters are different. */
0095         friend bool operator!=(const param_type& lhs, const param_type& rhs)
0096         {
0097             return !(lhs == rhs);
0098         }
0099     private:
0100         RealType _alpha;
0101         RealType _beta;
0102     };
0103
0104 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
0105     BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
0106 #endif
0107
0108     /**
0109      * Creates a new gamma_distribution with parameters "alpha" and "beta".
0110      *
0111      * Requires: alpha > 0 && beta > 0
0112      */
0113     explicit gamma_distribution(const result_type& alpha_arg = result_type(1.0),
0114                                 const result_type& beta_arg = result_type(1.0))
0115       : _exp(result_type(1)), _alpha(alpha_arg), _beta(beta_arg)
0116     {
0117         BOOST_ASSERT(_alpha > result_type(0));
0118         BOOST_ASSERT(_beta > result_type(0));
0119         init();
0120     }
0121
0122     /** Constructs a @c gamma_distribution from its parameters. */
0123     explicit gamma_distribution(const param_type& parm)
0124       : _exp(result_type(1)), _alpha(parm.alpha()), _beta(parm.beta())
0125     {
0126         init();
0127     }
0128
0129     // compiler-generated copy ctor and assignment operator are fine
0130
0131     /** Returns the "alpha" paramter of the distribution. */
0132     RealType alpha() const { return _alpha; }
0133     /** Returns the "beta" parameter of the distribution. */
0134     RealType beta() const { return _beta; }
0135     /** Returns the smallest value that the distribution can produce. */
0136     RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
0137     /* Returns the largest value that the distribution can produce. */
0138     RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
0139     { return (std::numeric_limits<RealType>::infinity)(); }
0140
0141     /** Returns the parameters of the distribution. */
0142     param_type param() const { return param_type(_alpha, _beta); }
0143     /** Sets the parameters of the distribution. */
0144     void param(const param_type& parm)
0145     {
0146         _alpha = parm.alpha();
0147         _beta = parm.beta();
0148         init();
0149     }
0150
0151     /**
0152      * Effects: Subsequent uses of the distribution do not depend
0153      * on values produced by any engine prior to invoking reset.
0154      */
0155     void reset() { _exp.reset(); }
0156
0157     /**
0158      * Returns a random variate distributed according to
0159      * the gamma distribution.
0160      */
0161     template<class Engine>
0162     result_type operator()(Engine& eng)
0163     {
0164 #ifndef BOOST_NO_STDC_NAMESPACE
0165         // allow for Koenig lookup
0166         using std::tan; using std::sqrt; using std::exp; using std::log;
0167         using std::pow;
0168 #endif
0169         if(_alpha == result_type(1)) {
0170             return _exp(eng) * _beta;
0171         } else if(_alpha > result_type(1)) {
0172             // Can we have a boost::mathconst please?
0173             const result_type pi = result_type(3.14159265358979323846);
0174             for(;;) {
0175                 result_type y = tan(pi * uniform_01<RealType>()(eng));
0176                 result_type x = sqrt(result_type(2)*_alpha-result_type(1))*y
0177                     + _alpha-result_type(1);
0178                 if(x <= result_type(0))
0179                     continue;
0180                 if(uniform_01<RealType>()(eng) >
0181                     (result_type(1)+y*y) * exp((_alpha-result_type(1))
0182                                                *log(x/(_alpha-result_type(1)))
0183                                                - sqrt(result_type(2)*_alpha
0184                                                       -result_type(1))*y))
0185                     continue;
0186                 return x * _beta;
0187             }
0188         } else /* alpha < 1.0 */ {
0189             for(;;) {
0190                 result_type u = uniform_01<RealType>()(eng);
0191                 result_type y = _exp(eng);
0192                 result_type x, q;
0193                 if(u < _p) {
0194                     x = exp(-y/_alpha);
0195                     q = _p*exp(-x);
0196                 } else {
0197                     x = result_type(1)+y;
0198                     q = _p + (result_type(1)-_p) * pow(x,_alpha-result_type(1));
0199                 }
0200                 if(u >= q)
0201                     continue;
0202                 return x * _beta;
0203             }
0204         }
0205     }
0206
0207     template<class URNG>
0208     RealType operator()(URNG& urng, const param_type& parm) const
0209     {
0210         return gamma_distribution(parm)(urng);
0211     }
0212
0213 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
0214     /** Writes a @c gamma_distribution to a @c std::ostream. */
0215     template<class CharT, class Traits>
0216     friend std::basic_ostream<CharT,Traits>&
0217     operator<<(std::basic_ostream<CharT,Traits>& os,
0218                const gamma_distribution& gd)
0219     {
0220         os << gd.param();
0221         return os;
0222     }
0223
0224     /** Reads a @c gamma_distribution from a @c std::istream. */
0225     template<class CharT, class Traits>
0226     friend std::basic_istream<CharT,Traits>&
0227     operator>>(std::basic_istream<CharT,Traits>& is, gamma_distribution& gd)
0228     {
0229         gd.read(is);
0230         return is;
0231     }
0232 #endif
0233
0234     /**
0235      * Returns true if the two distributions will produce identical
0236      * sequences of random variates given equal generators.
0237      */
0238     friend bool operator==(const gamma_distribution& lhs,
0239                            const gamma_distribution& rhs)
0240     {
0241         return lhs._alpha == rhs._alpha
0242             && lhs._beta == rhs._beta
0243             && lhs._exp == rhs._exp;
0244     }
0245
0246     /**
0247      * Returns true if the two distributions can produce different
0248      * sequences of random variates, given equal generators.
0249      */
0250     friend bool operator!=(const gamma_distribution& lhs,
0251                            const gamma_distribution& rhs)
0252     {
0253         return !(lhs == rhs);
0254     }
0255
0256 private:
0257     /// \cond hide_private_members
0258
0259     template<class CharT, class Traits>
0260     void read(std::basic_istream<CharT, Traits>& is)
0261     {
0262         param_type parm;
0263         if(is >> parm) {
0264             param(parm);
0265         }
0266     }
0267
0268     void init()
0269     {
0270 #ifndef BOOST_NO_STDC_NAMESPACE
0271         // allow for Koenig lookup
0272         using std::exp;
0273 #endif
0274         _p = exp(result_type(1)) / (_alpha + exp(result_type(1)));
0275     }
0276     /// \endcond
0277
0278     exponential_distribution<RealType> _exp;
0279     result_type _alpha;
0280     result_type _beta;
0281     // some data precomputed from the parameters
0282     result_type _p;
0283 };
0284
0285
0286 } // namespace random
0287
0288 using random::gamma_distribution;
0289
0290 } // namespace boost
0291
0292 #endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP