1 /* boost random/gamma_distribution.hpp header file
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)
9 * See http://www.boost.org for most recent version including documentation.
11 * $Id: gamma_distribution.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $
15 #ifndef BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
16 #define BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP
18 #include <boost/config/no_tr1/cmath.hpp>
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>
30 // The algorithm is taken from Knuth
33 * The gamma distribution is a continuous distribution with two
34 * parameters alpha and beta. It produces values > 0.
37 * \f$\displaystyle p(x) = x^{\alpha-1}\frac{e^{-x/\beta}}{\beta^\alpha\Gamma(\alpha)}\f$.
39 template<class RealType = double>
40 class gamma_distribution
43 typedef RealType input_type;
44 typedef RealType result_type;
49 typedef gamma_distribution distribution_type;
52 * Constructs a @c param_type object from the "alpha" and "beta"
55 * Requires: alpha > 0 && beta > 0
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)
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; }
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)
75 os << parm._alpha << ' ' << parm._beta;
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)
84 is >> parm._alpha >> std::ws >> parm._beta;
89 /** Returns true if the two sets of parameters are the same. */
90 friend bool operator==(const param_type& lhs, const param_type& rhs)
92 return lhs._alpha == rhs._alpha && lhs._beta == rhs._beta;
94 /** Returns true if the two sets fo parameters are different. */
95 friend bool operator!=(const param_type& lhs, const param_type& rhs)
104 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
105 BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer);
109 * Creates a new gamma_distribution with parameters "alpha" and "beta".
111 * Requires: alpha > 0 && beta > 0
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)
117 BOOST_ASSERT(_alpha > result_type(0));
118 BOOST_ASSERT(_beta > result_type(0));
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())
129 // compiler-generated copy ctor and assignment operator are fine
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)(); }
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)
146 _alpha = parm.alpha();
152 * Effects: Subsequent uses of the distribution do not depend
153 * on values produced by any engine prior to invoking reset.
155 void reset() { _exp.reset(); }
158 * Returns a random variate distributed according to
159 * the gamma distribution.
161 template<class Engine>
162 result_type operator()(Engine& eng)
164 #ifndef BOOST_NO_STDC_NAMESPACE
165 // allow for Koenig lookup
166 using std::tan; using std::sqrt; using std::exp; using std::log;
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);
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))
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
188 } else /* alpha < 1.0 */ {
190 result_type u = uniform_01<RealType>()(eng);
191 result_type y = _exp(eng);
197 x = result_type(1)+y;
198 q = _p + (result_type(1)-_p) * pow(x,_alpha-result_type(1));
208 RealType operator()(URNG& urng, const param_type& parm) const
210 return gamma_distribution(parm)(urng);
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)
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)
235 * Returns true if the two distributions will produce identical
236 * sequences of random variates given equal generators.
238 friend bool operator==(const gamma_distribution& lhs,
239 const gamma_distribution& rhs)
241 return lhs._alpha == rhs._alpha
242 && lhs._beta == rhs._beta
243 && lhs._exp == rhs._exp;
247 * Returns true if the two distributions can produce different
248 * sequences of random variates, given equal generators.
250 friend bool operator!=(const gamma_distribution& lhs,
251 const gamma_distribution& rhs)
253 return !(lhs == rhs);
257 /// \cond hide_private_members
259 template<class CharT, class Traits>
260 void read(std::basic_istream<CharT, Traits>& is)
270 #ifndef BOOST_NO_STDC_NAMESPACE
271 // allow for Koenig lookup
274 _p = exp(result_type(1)) / (_alpha + exp(result_type(1)));
278 exponential_distribution<RealType> _exp;
281 // some data precomputed from the parameters
286 } // namespace random
288 using random::gamma_distribution;
292 #endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP