diff options
Diffstat (limited to 'src/third_party/boost-1.69.0/boost/random/gamma_distribution.hpp')
| -rw-r--r-- | src/third_party/boost-1.69.0/boost/random/gamma_distribution.hpp | 292 |
1 files changed, 292 insertions, 0 deletions
diff --git a/src/third_party/boost-1.69.0/boost/random/gamma_distribution.hpp b/src/third_party/boost-1.69.0/boost/random/gamma_distribution.hpp new file mode 100644 index 00000000000..c9bda4ca9d3 --- /dev/null +++ b/src/third_party/boost-1.69.0/boost/random/gamma_distribution.hpp @@ -0,0 +1,292 @@ +/* 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 BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP +#define BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP + +#include <boost/config/no_tr1/cmath.hpp> +#include <istream> +#include <iosfwd> +#include <boost/assert.hpp> +#include <boost/limits.hpp> +#include <boost/static_assert.hpp> +#include <boost/random/detail/config.hpp> +#include <boost/random/exponential_distribution.hpp> + +namespace boost { +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 RealType = double> +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<class CharT, class Traits> + friend std::basic_ostream<CharT, Traits>& + operator<<(std::basic_ostream<CharT, Traits>& os, + const param_type& parm) + { + os << parm._alpha << ' ' << parm._beta; + return os; + } + + /** Reads the parameters from a @c std::istream. */ + template<class CharT, class Traits> + friend std::basic_istream<CharT, Traits>& + operator>>(std::basic_istream<CharT, Traits>& 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<RealType>::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<RealType>::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<class Engine> + 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 * uniform_01<RealType>()(eng)); + result_type x = sqrt(result_type(2)*_alpha-result_type(1))*y + + _alpha-result_type(1); + if(x <= result_type(0)) + continue; + if(uniform_01<RealType>()(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 = uniform_01<RealType>()(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<class URNG> + 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<class CharT, class Traits> + friend std::basic_ostream<CharT,Traits>& + operator<<(std::basic_ostream<CharT,Traits>& os, + const gamma_distribution& gd) + { + os << gd.param(); + return os; + } + + /** Reads a @c gamma_distribution from a @c std::istream. */ + template<class CharT, class Traits> + friend std::basic_istream<CharT,Traits>& + operator>>(std::basic_istream<CharT,Traits>& 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<class CharT, class Traits> + void read(std::basic_istream<CharT, Traits>& 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 + + exponential_distribution<RealType> _exp; + result_type _alpha; + result_type _beta; + // some data precomputed from the parameters + result_type _p; +}; + + +} // namespace random + +using random::gamma_distribution; + +} // namespace boost + +#endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP |