Commit 7ec3e39b by Manuel Günther

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