Chris@16: /* boost random/gamma_distribution.hpp header file
Chris@16:  *
Chris@16:  * Copyright Jens Maurer 2002
Chris@16:  * Copyright Steven Watanabe 2010
Chris@16:  * Distributed under the Boost Software License, Version 1.0. (See
Chris@16:  * accompanying file LICENSE_1_0.txt or copy at
Chris@16:  * http://www.boost.org/LICENSE_1_0.txt)
Chris@16:  *
Chris@16:  * See http://www.boost.org for most recent version including documentation.
Chris@16:  *
Chris@101:  * $Id$
Chris@16:  *
Chris@16:  */
Chris@16: 
Chris@16: #ifndef BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
Chris@16: #define BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
Chris@16: 
Chris@16: #include <boost/config/no_tr1/cmath.hpp>
Chris@16: #include <istream>
Chris@16: #include <iosfwd>
Chris@16: #include <boost/assert.hpp>
Chris@16: #include <boost/limits.hpp>
Chris@16: #include <boost/static_assert.hpp>
Chris@16: #include <boost/random/detail/config.hpp>
Chris@16: #include <boost/random/exponential_distribution.hpp>
Chris@16: 
Chris@16: namespace boost {
Chris@16: namespace random {
Chris@16: 
Chris@16: // The algorithm is taken from Knuth
Chris@16: 
Chris@16: /**
Chris@16:  * The gamma distribution is a continuous distribution with two
Chris@16:  * parameters alpha and beta.  It produces values > 0.
Chris@16:  *
Chris@16:  * It has
Chris@16:  * \f$\displaystyle p(x) = x^{\alpha-1}\frac{e^{-x/\beta}}{\beta^\alpha\Gamma(\alpha)}\f$.
Chris@16:  */
Chris@16: template<class RealType = double>
Chris@16: class gamma_distribution
Chris@16: {
Chris@16: public:
Chris@16:     typedef RealType input_type;
Chris@16:     typedef RealType result_type;
Chris@16: 
Chris@16:     class param_type
Chris@16:     {
Chris@16:     public:
Chris@16:         typedef gamma_distribution distribution_type;
Chris@16: 
Chris@16:         /**
Chris@16:          * Constructs a @c param_type object from the "alpha" and "beta"
Chris@16:          * parameters.
Chris@16:          *
Chris@16:          * Requires: alpha > 0 && beta > 0
Chris@16:          */
Chris@16:         param_type(const RealType& alpha_arg = RealType(1.0),
Chris@16:                    const RealType& beta_arg = RealType(1.0))
Chris@16:           : _alpha(alpha_arg), _beta(beta_arg)
Chris@16:         {
Chris@16:         }
Chris@16: 
Chris@16:         /** Returns the "alpha" parameter of the distribution. */
Chris@16:         RealType alpha() const { return _alpha; }
Chris@16:         /** Returns the "beta" parameter of the distribution. */
Chris@16:         RealType beta() const { return _beta; }
Chris@16: 
Chris@16: #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
Chris@16:         /** Writes the parameters to a @c std::ostream. */
Chris@16:         template<class CharT, class Traits>
Chris@16:         friend std::basic_ostream<CharT, Traits>&
Chris@16:         operator<<(std::basic_ostream<CharT, Traits>& os,
Chris@16:                    const param_type& parm)
Chris@16:         {
Chris@16:             os << parm._alpha << ' ' << parm._beta;
Chris@16:             return os;
Chris@16:         }
Chris@16:         
Chris@16:         /** Reads the parameters from a @c std::istream. */
Chris@16:         template<class CharT, class Traits>
Chris@16:         friend std::basic_istream<CharT, Traits>&
Chris@16:         operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm)
Chris@16:         {
Chris@16:             is >> parm._alpha >> std::ws >> parm._beta;
Chris@16:             return is;
Chris@16:         }
Chris@16: #endif
Chris@16: 
Chris@16:         /** Returns true if the two sets of parameters are the same. */
Chris@16:         friend bool operator==(const param_type& lhs, const param_type& rhs)
Chris@16:         {
Chris@16:             return lhs._alpha == rhs._alpha && lhs._beta == rhs._beta;
Chris@16:         }
Chris@16:         /** Returns true if the two sets fo parameters are different. */
Chris@16:         friend bool operator!=(const param_type& lhs, const param_type& rhs)
Chris@16:         {
Chris@16:             return !(lhs == rhs);
Chris@16:         }
Chris@16:     private:
Chris@16:         RealType _alpha;
Chris@16:         RealType _beta;
Chris@16:     };
Chris@16: 
Chris@16: #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
Chris@16:     BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
Chris@16: #endif
Chris@16: 
Chris@16:     /**
Chris@16:      * Creates a new gamma_distribution with parameters "alpha" and "beta".
Chris@16:      *
Chris@16:      * Requires: alpha > 0 && beta > 0
Chris@16:      */
Chris@16:     explicit gamma_distribution(const result_type& alpha_arg = result_type(1.0),
Chris@16:                                 const result_type& beta_arg = result_type(1.0))
Chris@16:       : _exp(result_type(1)), _alpha(alpha_arg), _beta(beta_arg)
Chris@16:     {
Chris@16:         BOOST_ASSERT(_alpha > result_type(0));
Chris@16:         BOOST_ASSERT(_beta > result_type(0));
Chris@16:         init();
Chris@16:     }
Chris@16: 
Chris@16:     /** Constructs a @c gamma_distribution from its parameters. */
Chris@16:     explicit gamma_distribution(const param_type& parm)
Chris@16:       : _exp(result_type(1)), _alpha(parm.alpha()), _beta(parm.beta())
Chris@16:     {
Chris@16:         init();
Chris@16:     }
Chris@16: 
Chris@16:     // compiler-generated copy ctor and assignment operator are fine
Chris@16: 
Chris@16:     /** Returns the "alpha" paramter of the distribution. */
Chris@16:     RealType alpha() const { return _alpha; }
Chris@16:     /** Returns the "beta" parameter of the distribution. */
Chris@16:     RealType beta() const { return _beta; }
Chris@16:     /** Returns the smallest value that the distribution can produce. */
Chris@16:     RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }
Chris@16:     /* Returns the largest value that the distribution can produce. */
Chris@16:     RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const
Chris@16:     { return (std::numeric_limits<RealType>::infinity)(); }
Chris@16: 
Chris@16:     /** Returns the parameters of the distribution. */
Chris@16:     param_type param() const { return param_type(_alpha, _beta); }
Chris@16:     /** Sets the parameters of the distribution. */
Chris@16:     void param(const param_type& parm)
Chris@16:     {
Chris@16:         _alpha = parm.alpha();
Chris@16:         _beta = parm.beta();
Chris@16:         init();
Chris@16:     }
Chris@16:     
Chris@16:     /**
Chris@16:      * Effects: Subsequent uses of the distribution do not depend
Chris@16:      * on values produced by any engine prior to invoking reset.
Chris@16:      */
Chris@16:     void reset() { _exp.reset(); }
Chris@16: 
Chris@16:     /**
Chris@16:      * Returns a random variate distributed according to
Chris@16:      * the gamma distribution.
Chris@16:      */
Chris@16:     template<class Engine>
Chris@16:     result_type operator()(Engine& eng)
Chris@16:     {
Chris@16: #ifndef BOOST_NO_STDC_NAMESPACE
Chris@16:         // allow for Koenig lookup
Chris@16:         using std::tan; using std::sqrt; using std::exp; using std::log;
Chris@16:         using std::pow;
Chris@16: #endif
Chris@16:         if(_alpha == result_type(1)) {
Chris@16:             return _exp(eng) * _beta;
Chris@16:         } else if(_alpha > result_type(1)) {
Chris@16:             // Can we have a boost::mathconst please?
Chris@16:             const result_type pi = result_type(3.14159265358979323846);
Chris@16:             for(;;) {
Chris@16:                 result_type y = tan(pi * uniform_01<RealType>()(eng));
Chris@16:                 result_type x = sqrt(result_type(2)*_alpha-result_type(1))*y
Chris@16:                     + _alpha-result_type(1);
Chris@16:                 if(x <= result_type(0))
Chris@16:                     continue;
Chris@16:                 if(uniform_01<RealType>()(eng) >
Chris@16:                     (result_type(1)+y*y) * exp((_alpha-result_type(1))
Chris@16:                                                *log(x/(_alpha-result_type(1)))
Chris@16:                                                - sqrt(result_type(2)*_alpha
Chris@16:                                                       -result_type(1))*y))
Chris@16:                     continue;
Chris@16:                 return x * _beta;
Chris@16:             }
Chris@16:         } else /* alpha < 1.0 */ {
Chris@16:             for(;;) {
Chris@16:                 result_type u = uniform_01<RealType>()(eng);
Chris@16:                 result_type y = _exp(eng);
Chris@16:                 result_type x, q;
Chris@16:                 if(u < _p) {
Chris@16:                     x = exp(-y/_alpha);
Chris@16:                     q = _p*exp(-x);
Chris@16:                 } else {
Chris@16:                     x = result_type(1)+y;
Chris@16:                     q = _p + (result_type(1)-_p) * pow(x,_alpha-result_type(1));
Chris@16:                 }
Chris@16:                 if(u >= q)
Chris@16:                     continue;
Chris@16:                 return x * _beta;
Chris@16:             }
Chris@16:         }
Chris@16:     }
Chris@16: 
Chris@16:     template<class URNG>
Chris@16:     RealType operator()(URNG& urng, const param_type& parm) const
Chris@16:     {
Chris@16:         return gamma_distribution(parm)(urng);
Chris@16:     }
Chris@16: 
Chris@16: #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS
Chris@16:     /** Writes a @c gamma_distribution to a @c std::ostream. */
Chris@16:     template<class CharT, class Traits>
Chris@16:     friend std::basic_ostream<CharT,Traits>&
Chris@16:     operator<<(std::basic_ostream<CharT,Traits>& os,
Chris@16:                const gamma_distribution& gd)
Chris@16:     {
Chris@16:         os << gd.param();
Chris@16:         return os;
Chris@16:     }
Chris@16:     
Chris@16:     /** Reads a @c gamma_distribution from a @c std::istream. */
Chris@16:     template<class CharT, class Traits>
Chris@16:     friend std::basic_istream<CharT,Traits>&
Chris@16:     operator>>(std::basic_istream<CharT,Traits>& is, gamma_distribution& gd)
Chris@16:     {
Chris@16:         gd.read(is);
Chris@16:         return is;
Chris@16:     }
Chris@16: #endif
Chris@16: 
Chris@16:     /**
Chris@16:      * Returns true if the two distributions will produce identical
Chris@16:      * sequences of random variates given equal generators.
Chris@16:      */
Chris@16:     friend bool operator==(const gamma_distribution& lhs,
Chris@16:                            const gamma_distribution& rhs)
Chris@16:     {
Chris@16:         return lhs._alpha == rhs._alpha
Chris@16:             && lhs._beta == rhs._beta
Chris@16:             && lhs._exp == rhs._exp;
Chris@16:     }
Chris@16: 
Chris@16:     /**
Chris@16:      * Returns true if the two distributions can produce different
Chris@16:      * sequences of random variates, given equal generators.
Chris@16:      */
Chris@16:     friend bool operator!=(const gamma_distribution& lhs,
Chris@16:                            const gamma_distribution& rhs)
Chris@16:     {
Chris@16:         return !(lhs == rhs);
Chris@16:     }
Chris@16: 
Chris@16: private:
Chris@16:     /// \cond hide_private_members
Chris@16: 
Chris@16:     template<class CharT, class Traits>
Chris@16:     void read(std::basic_istream<CharT, Traits>& is)
Chris@16:     {
Chris@16:         param_type parm;
Chris@16:         if(is >> parm) {
Chris@16:             param(parm);
Chris@16:         }
Chris@16:     }
Chris@16: 
Chris@16:     void init()
Chris@16:     {
Chris@16: #ifndef BOOST_NO_STDC_NAMESPACE
Chris@16:         // allow for Koenig lookup
Chris@16:         using std::exp;
Chris@16: #endif
Chris@16:         _p = exp(result_type(1)) / (_alpha + exp(result_type(1)));
Chris@16:     }
Chris@16:     /// \endcond
Chris@16: 
Chris@16:     exponential_distribution<RealType> _exp;
Chris@16:     result_type _alpha;
Chris@16:     result_type _beta;
Chris@16:     // some data precomputed from the parameters
Chris@16:     result_type _p;
Chris@16: };
Chris@16: 
Chris@16: 
Chris@16: } // namespace random
Chris@16: 
Chris@16: using random::gamma_distribution;
Chris@16: 
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP