1 /* boost random/gamma_distribution.hpp header file 2 * 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) 8 * 9 * See http://www.boost.org for most recent version including documentation. 10 * 11 * $Id$ 12 * 13 */ 14 15 #ifndef BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP 16 #define BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP 17 18 #include <boost/config/no_tr1/cmath.hpp> 19 #include <istream> 20 #include <iosfwd> 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> 26 27 namespace boost { 28 namespace random { 29 30 // The algorithm is taken from Knuth 31 32 /** 33 * The gamma distribution is a continuous distribution with two 34 * parameters alpha and beta. It produces values > 0. 35 * 36 * It has 37 * \f$\displaystyle p(x) = x^{\alpha-1}\frac{e^{-x/\beta}}{\beta^\alpha\Gamma(\alpha)}\f$. 38 */ 39 template<class RealType = double> 40 class gamma_distribution 41 { 42 public: 43 typedef RealType input_type; 44 typedef RealType result_type; 45 46 class param_type 47 { 48 public: 49 typedef gamma_distribution distribution_type; 50 51 /** 52 * Constructs a @c param_type object from the "alpha" and "beta" 53 * parameters. 54 * 55 * Requires: alpha > 0 && beta > 0 56 */ param_type(const RealType & alpha_arg=RealType (1.0),const RealType & beta_arg=RealType (1.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) 60 { 61 } 62 63 /** Returns the "alpha" parameter of the distribution. */ alpha() const64 RealType alpha() const { return _alpha; } 65 /** Returns the "beta" parameter of the distribution. */ beta() const66 RealType beta() const { return _beta; } 67 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>& operator <<(std::basic_ostream<CharT,Traits> & os,const param_type & parm)72 operator<<(std::basic_ostream<CharT, Traits>& os, 73 const param_type& parm) 74 { 75 os << parm._alpha << ' ' << parm._beta; 76 return os; 77 } 78 79 /** Reads the parameters from a @c std::istream. */ 80 template<class CharT, class Traits> 81 friend std::basic_istream<CharT, Traits>& operator >>(std::basic_istream<CharT,Traits> & is,param_type & parm)82 operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm) 83 { 84 is >> parm._alpha >> std::ws >> parm._beta; 85 return is; 86 } 87 #endif 88 89 /** Returns true if the two sets of parameters are the same. */ operator ==(const param_type & lhs,const param_type & rhs)90 friend bool operator==(const param_type& lhs, const param_type& rhs) 91 { 92 return lhs._alpha == rhs._alpha && lhs._beta == rhs._beta; 93 } 94 /** Returns true if the two sets fo parameters are different. */ operator !=(const param_type & lhs,const param_type & rhs)95 friend bool operator!=(const param_type& lhs, const param_type& rhs) 96 { 97 return !(lhs == rhs); 98 } 99 private: 100 RealType _alpha; 101 RealType _beta; 102 }; 103 104 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS 105 BOOST_STATIC_ASSERT(!std::numeric_limits<RealType>::is_integer); 106 #endif 107 108 /** 109 * Creates a new gamma_distribution with parameters "alpha" and "beta". 110 * 111 * Requires: alpha > 0 && beta > 0 112 */ gamma_distribution(const result_type & alpha_arg=result_type (1.0),const result_type & beta_arg=result_type (1.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) 116 { 117 BOOST_ASSERT(_alpha > result_type(0)); 118 BOOST_ASSERT(_beta > result_type(0)); 119 init(); 120 } 121 122 /** Constructs a @c gamma_distribution from its parameters. */ gamma_distribution(const param_type & parm)123 explicit gamma_distribution(const param_type& parm) 124 : _exp(result_type(1)), _alpha(parm.alpha()), _beta(parm.beta()) 125 { 126 init(); 127 } 128 129 // compiler-generated copy ctor and assignment operator are fine 130 131 /** Returns the "alpha" paramter of the distribution. */ alpha() const132 RealType alpha() const { return _alpha; } 133 /** Returns the "beta" parameter of the distribution. */ beta() const134 RealType beta() const { return _beta; } 135 /** Returns the smallest value that the distribution can produce. */ BOOST_PREVENT_MACRO_SUBSTITUTION() const136 RealType min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; } 137 /* Returns the largest value that the distribution can produce. */ BOOST_PREVENT_MACRO_SUBSTITUTION() const138 RealType max BOOST_PREVENT_MACRO_SUBSTITUTION () const 139 { return (std::numeric_limits<RealType>::infinity)(); } 140 141 /** Returns the parameters of the distribution. */ param() const142 param_type param() const { return param_type(_alpha, _beta); } 143 /** Sets the parameters of the distribution. */ param(const param_type & parm)144 void param(const param_type& parm) 145 { 146 _alpha = parm.alpha(); 147 _beta = parm.beta(); 148 init(); 149 } 150 151 /** 152 * Effects: Subsequent uses of the distribution do not depend 153 * on values produced by any engine prior to invoking reset. 154 */ reset()155 void reset() { _exp.reset(); } 156 157 /** 158 * Returns a random variate distributed according to 159 * the gamma distribution. 160 */ 161 template<class Engine> operator ()(Engine & eng)162 result_type operator()(Engine& eng) 163 { 164 #ifndef BOOST_NO_STDC_NAMESPACE 165 // allow for Koenig lookup 166 using std::tan; using std::sqrt; using std::exp; using std::log; 167 using std::pow; 168 #endif 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); 174 for(;;) { 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)) 179 continue; 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 184 -result_type(1))*y)) 185 continue; 186 return x * _beta; 187 } 188 } else /* alpha < 1.0 */ { 189 for(;;) { 190 result_type u = uniform_01<RealType>()(eng); 191 result_type y = _exp(eng); 192 result_type x, q; 193 if(u < _p) { 194 x = exp(-y/_alpha); 195 q = _p*exp(-x); 196 } else { 197 x = result_type(1)+y; 198 q = _p + (result_type(1)-_p) * pow(x,_alpha-result_type(1)); 199 } 200 if(u >= q) 201 continue; 202 return x * _beta; 203 } 204 } 205 } 206 207 template<class URNG> operator ()(URNG & urng,const param_type & parm) const208 RealType operator()(URNG& urng, const param_type& parm) const 209 { 210 return gamma_distribution(parm)(urng); 211 } 212 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>& operator <<(std::basic_ostream<CharT,Traits> & os,const gamma_distribution & gd)217 operator<<(std::basic_ostream<CharT,Traits>& os, 218 const gamma_distribution& gd) 219 { 220 os << gd.param(); 221 return os; 222 } 223 224 /** Reads a @c gamma_distribution from a @c std::istream. */ 225 template<class CharT, class Traits> 226 friend std::basic_istream<CharT,Traits>& operator >>(std::basic_istream<CharT,Traits> & is,gamma_distribution & gd)227 operator>>(std::basic_istream<CharT,Traits>& is, gamma_distribution& gd) 228 { 229 gd.read(is); 230 return is; 231 } 232 #endif 233 234 /** 235 * Returns true if the two distributions will produce identical 236 * sequences of random variates given equal generators. 237 */ operator ==(const gamma_distribution & lhs,const gamma_distribution & rhs)238 friend bool operator==(const gamma_distribution& lhs, 239 const gamma_distribution& rhs) 240 { 241 return lhs._alpha == rhs._alpha 242 && lhs._beta == rhs._beta 243 && lhs._exp == rhs._exp; 244 } 245 246 /** 247 * Returns true if the two distributions can produce different 248 * sequences of random variates, given equal generators. 249 */ operator !=(const gamma_distribution & lhs,const gamma_distribution & rhs)250 friend bool operator!=(const gamma_distribution& lhs, 251 const gamma_distribution& rhs) 252 { 253 return !(lhs == rhs); 254 } 255 256 private: 257 /// \cond hide_private_members 258 259 template<class CharT, class Traits> read(std::basic_istream<CharT,Traits> & is)260 void read(std::basic_istream<CharT, Traits>& is) 261 { 262 param_type parm; 263 if(is >> parm) { 264 param(parm); 265 } 266 } 267 init()268 void init() 269 { 270 #ifndef BOOST_NO_STDC_NAMESPACE 271 // allow for Koenig lookup 272 using std::exp; 273 #endif 274 _p = exp(result_type(1)) / (_alpha + exp(result_type(1))); 275 } 276 /// \endcond 277 278 exponential_distribution<RealType> _exp; 279 result_type _alpha; 280 result_type _beta; 281 // some data precomputed from the parameters 282 result_type _p; 283 }; 284 285 286 } // namespace random 287 288 using random::gamma_distribution; 289 290 } // namespace boost 291 292 #endif // BOOST_RANDOM_GAMMA_DISTRIBUTION_HPP 293