Ticket #6135: gamma_distribution.patch
File gamma_distribution.patch, 5.3 KB (added by , 11 years ago) |
---|
-
.hpp
old new 22 22 #include <boost/limits.hpp> 23 23 #include <boost/static_assert.hpp> 24 24 #include <boost/random/detail/config.hpp> 25 #include <boost/random/ exponential_distribution.hpp>25 #include <boost/random/normal_distribution.hpp> 26 26 27 27 namespace boost { 28 28 namespace random { 29 29 30 // The algorithm is taken from Knuth 30 // The algorithm is taken from Marsaglia and Tsang, "A Simple Method 31 // for generating gamma variables", ACM Transactions on Mathematical 32 // Software, Vol 26, No 3 (2000), p363-372, as implemented in the 33 // GNU Scientific Library, v1.15 31 34 32 35 /** 33 36 * The gamma distribution is a continuous distribution with two … … 112 115 */ 113 116 explicit gamma_distribution(const result_type& alpha_arg = result_type(1.0), 114 117 const result_type& beta_arg = result_type(1.0)) 115 : _ exp(result_type(1)), _alpha(alpha_arg), _beta(beta_arg)118 : _normal(), _alpha(alpha_arg), _beta(beta_arg) 116 119 { 117 120 BOOST_ASSERT(_alpha > result_type(0)); 118 121 BOOST_ASSERT(_beta > result_type(0)); … … 121 124 122 125 /** Constructs a @c gamma_distribution from its parameters. */ 123 126 explicit gamma_distribution(const param_type& parm) 124 : _ exp(result_type(1)), _alpha(parm.alpha()), _beta(parm.beta())127 : _normal(), _alpha(parm.alpha()), _beta(parm.beta()) 125 128 { 126 129 init(); 127 130 } … … 152 155 * Effects: Subsequent uses of the distribution do not depend 153 156 * on values produced by any engine prior to invoking reset. 154 157 */ 155 void reset() { _ exp.reset(); }158 void reset() { _normal.reset(); } 156 159 157 160 /** 158 161 * Returns a random variate distributed according to … … 163 166 { 164 167 #ifndef BOOST_NO_STDC_NAMESPACE 165 168 // allow for Koenig lookup 166 using std::tan; using std::sqrt; using std::exp; using std::log; 167 using std::pow; 169 using std::log; using std::pow; 168 170 #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 } 171 if(_alpha < result_type(1)) { 172 result_type u; 173 do u = uniform_01<RealType>()(eng); 174 while(u <= result_type(0)); 175 return gamma_distribution<RealType>(_alpha + result_type(1), _beta)(eng) * 176 pow(u, result_type(1) / _alpha); 204 177 } 178 179 result_type v; 180 for(;;) { 181 result_type x; 182 do { 183 x = _normal(eng); 184 v = result_type(1) + _c * x; 185 } while(v <= result_type(0)); 186 187 v = v * v * v; 188 result_type u; 189 do u = uniform_01<RealType>()(eng); 190 while(u <= result_type(0)); 191 192 if (u < result_type(1) - result_type(0.0331) * x * x * x * x) 193 break; 194 195 if (log(u) < result_type(0.5) * x * x + _d * (result_type(1) - v + log(v))) 196 break; 197 } 198 return _beta * _d * v; 205 199 } 206 200 207 201 template<class URNG> … … 240 234 { 241 235 return lhs._alpha == rhs._alpha 242 236 && lhs._beta == rhs._beta 243 && lhs._ exp == rhs._exp;237 && lhs._normal == rhs._normal; 244 238 } 245 239 246 240 /** … … 269 263 { 270 264 #ifndef BOOST_NO_STDC_NAMESPACE 271 265 // allow for Koenig lookup 272 using std:: exp;266 using std::sqrt; 273 267 #endif 274 _p = exp(result_type(1)) / (_alpha + exp(result_type(1))); 268 _d = _alpha - result_type(1) / result_type(3); 269 _c = (result_type(1) / result_type(3)) / sqrt(_d); 275 270 } 276 271 /// \endcond 277 272 278 exponential_distribution<RealType> _exp;273 normal_distribution<RealType> _normal; 279 274 result_type _alpha; 280 275 result_type _beta; 281 276 // some data precomputed from the parameters 282 result_type _p; 277 result_type _d; 278 result_type _c; 283 279 }; 284 280 285 281