1 /* boost random/poisson_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_POISSON_DISTRIBUTION_HPP 16 #define BOOST_RANDOM_POISSON_DISTRIBUTION_HPP 17 18 #include <boost/config/no_tr1/cmath.hpp> 19 #include <cstdlib> 20 #include <iosfwd> 21 #include <boost/assert.hpp> 22 #include <boost/limits.hpp> 23 #include <boost/random/uniform_01.hpp> 24 #include <boost/random/detail/config.hpp> 25 26 #include <boost/random/detail/disable_warnings.hpp> 27 28 namespace boost { 29 namespace random { 30 31 namespace detail { 32 33 template<class RealType> 34 struct poisson_table { 35 static RealType value[10]; 36 }; 37 38 template<class RealType> 39 RealType poisson_table<RealType>::value[10] = { 40 0.0, 41 0.0, 42 0.69314718055994529, 43 1.7917594692280550, 44 3.1780538303479458, 45 4.7874917427820458, 46 6.5792512120101012, 47 8.5251613610654147, 48 10.604602902745251, 49 12.801827480081469 50 }; 51 52 } 53 54 /** 55 * An instantiation of the class template @c poisson_distribution is a 56 * model of \random_distribution. The poisson distribution has 57 * \f$p(i) = \frac{e^{-\lambda}\lambda^i}{i!}\f$ 58 * 59 * This implementation is based on the PTRD algorithm described 60 * 61 * @blockquote 62 * "The transformed rejection method for generating Poisson random variables", 63 * Wolfgang Hormann, Insurance: Mathematics and Economics 64 * Volume 12, Issue 1, February 1993, Pages 39-45 65 * @endblockquote 66 */ 67 template<class IntType = int, class RealType = double> 68 class poisson_distribution { 69 public: 70 typedef IntType result_type; 71 typedef RealType input_type; 72 73 class param_type { 74 public: 75 typedef poisson_distribution distribution_type; 76 /** 77 * Construct a param_type object with the parameter "mean" 78 * 79 * Requires: mean > 0 80 */ param_type(RealType mean_arg=RealType (1))81 explicit param_type(RealType mean_arg = RealType(1)) 82 : _mean(mean_arg) 83 { 84 BOOST_ASSERT(_mean > 0); 85 } 86 /* Returns the "mean" parameter of the distribution. */ mean() const87 RealType mean() const { return _mean; } 88 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS 89 /** Writes the parameters of the distribution to a @c std::ostream. */ 90 template<class CharT, class Traits> 91 friend std::basic_ostream<CharT, Traits>& operator <<(std::basic_ostream<CharT,Traits> & os,const param_type & parm)92 operator<<(std::basic_ostream<CharT, Traits>& os, 93 const param_type& parm) 94 { 95 os << parm._mean; 96 return os; 97 } 98 99 /** Reads the parameters of the distribution from a @c std::istream. */ 100 template<class CharT, class Traits> 101 friend std::basic_istream<CharT, Traits>& operator >>(std::basic_istream<CharT,Traits> & is,param_type & parm)102 operator>>(std::basic_istream<CharT, Traits>& is, param_type& parm) 103 { 104 is >> parm._mean; 105 return is; 106 } 107 #endif 108 /** Returns true if the parameters have the same values. */ operator ==(const param_type & lhs,const param_type & rhs)109 friend bool operator==(const param_type& lhs, const param_type& rhs) 110 { 111 return lhs._mean == rhs._mean; 112 } 113 /** Returns true if the parameters have different values. */ operator !=(const param_type & lhs,const param_type & rhs)114 friend bool operator!=(const param_type& lhs, const param_type& rhs) 115 { 116 return !(lhs == rhs); 117 } 118 private: 119 RealType _mean; 120 }; 121 122 /** 123 * Constructs a @c poisson_distribution with the parameter @c mean. 124 * 125 * Requires: mean > 0 126 */ poisson_distribution(RealType mean_arg=RealType (1))127 explicit poisson_distribution(RealType mean_arg = RealType(1)) 128 : _mean(mean_arg) 129 { 130 BOOST_ASSERT(_mean > 0); 131 init(); 132 } 133 134 /** 135 * Construct an @c poisson_distribution object from the 136 * parameters. 137 */ poisson_distribution(const param_type & parm)138 explicit poisson_distribution(const param_type& parm) 139 : _mean(parm.mean()) 140 { 141 init(); 142 } 143 144 /** 145 * Returns a random variate distributed according to the 146 * poisson distribution. 147 */ 148 template<class URNG> operator ()(URNG & urng) const149 IntType operator()(URNG& urng) const 150 { 151 if(use_inversion()) { 152 return invert(urng); 153 } else { 154 return generate(urng); 155 } 156 } 157 158 /** 159 * Returns a random variate distributed according to the 160 * poisson distribution with parameters specified by param. 161 */ 162 template<class URNG> operator ()(URNG & urng,const param_type & parm) const163 IntType operator()(URNG& urng, const param_type& parm) const 164 { 165 return poisson_distribution(parm)(urng); 166 } 167 168 /** Returns the "mean" parameter of the distribution. */ mean() const169 RealType mean() const { return _mean; } 170 171 /** Returns the smallest value that the distribution can produce. */ BOOST_PREVENT_MACRO_SUBSTITUTION() const172 IntType min BOOST_PREVENT_MACRO_SUBSTITUTION() const { return 0; } 173 /** Returns the largest value that the distribution can produce. */ BOOST_PREVENT_MACRO_SUBSTITUTION() const174 IntType max BOOST_PREVENT_MACRO_SUBSTITUTION() const 175 { return (std::numeric_limits<IntType>::max)(); } 176 177 /** Returns the parameters of the distribution. */ param() const178 param_type param() const { return param_type(_mean); } 179 /** Sets parameters of the distribution. */ param(const param_type & parm)180 void param(const param_type& parm) 181 { 182 _mean = parm.mean(); 183 init(); 184 } 185 186 /** 187 * Effects: Subsequent uses of the distribution do not depend 188 * on values produced by any engine prior to invoking reset. 189 */ reset()190 void reset() { } 191 192 #ifndef BOOST_RANDOM_NO_STREAM_OPERATORS 193 /** Writes the parameters of the distribution to a @c std::ostream. */ 194 template<class CharT, class Traits> 195 friend std::basic_ostream<CharT,Traits>& operator <<(std::basic_ostream<CharT,Traits> & os,const poisson_distribution & pd)196 operator<<(std::basic_ostream<CharT,Traits>& os, 197 const poisson_distribution& pd) 198 { 199 os << pd.param(); 200 return os; 201 } 202 203 /** Reads the parameters of the distribution from a @c std::istream. */ 204 template<class CharT, class Traits> 205 friend std::basic_istream<CharT,Traits>& operator >>(std::basic_istream<CharT,Traits> & is,poisson_distribution & pd)206 operator>>(std::basic_istream<CharT,Traits>& is, poisson_distribution& pd) 207 { 208 pd.read(is); 209 return is; 210 } 211 #endif 212 213 /** Returns true if the two distributions will produce the same 214 sequence of values, given equal generators. */ operator ==(const poisson_distribution & lhs,const poisson_distribution & rhs)215 friend bool operator==(const poisson_distribution& lhs, 216 const poisson_distribution& rhs) 217 { 218 return lhs._mean == rhs._mean; 219 } 220 /** Returns true if the two distributions could produce different 221 sequences of values, given equal generators. */ operator !=(const poisson_distribution & lhs,const poisson_distribution & rhs)222 friend bool operator!=(const poisson_distribution& lhs, 223 const poisson_distribution& rhs) 224 { 225 return !(lhs == rhs); 226 } 227 228 private: 229 230 /// @cond show_private 231 232 template<class CharT, class Traits> read(std::basic_istream<CharT,Traits> & is)233 void read(std::basic_istream<CharT, Traits>& is) { 234 param_type parm; 235 if(is >> parm) { 236 param(parm); 237 } 238 } 239 use_inversion() const240 bool use_inversion() const 241 { 242 return _mean < 10; 243 } 244 log_factorial(IntType k)245 static RealType log_factorial(IntType k) 246 { 247 BOOST_ASSERT(k >= 0); 248 BOOST_ASSERT(k < 10); 249 return detail::poisson_table<RealType>::value[k]; 250 } 251 init()252 void init() 253 { 254 using std::sqrt; 255 using std::exp; 256 257 if(use_inversion()) { 258 _u._exp_mean = exp(-_mean); 259 } else { 260 _u._ptrd.smu = sqrt(_mean); 261 _u._ptrd.b = 0.931 + 2.53 * _u._ptrd.smu; 262 _u._ptrd.a = -0.059 + 0.02483 * _u._ptrd.b; 263 _u._ptrd.inv_alpha = 1.1239 + 1.1328 / (_u._ptrd.b - 3.4); 264 _u._ptrd.v_r = 0.9277 - 3.6224 / (_u._ptrd.b - 2); 265 } 266 } 267 268 template<class URNG> generate(URNG & urng) const269 IntType generate(URNG& urng) const 270 { 271 using std::floor; 272 using std::abs; 273 using std::log; 274 275 while(true) { 276 RealType u; 277 RealType v = uniform_01<RealType>()(urng); 278 if(v <= 0.86 * _u._ptrd.v_r) { 279 u = v / _u._ptrd.v_r - 0.43; 280 return static_cast<IntType>(floor( 281 (2*_u._ptrd.a/(0.5-abs(u)) + _u._ptrd.b)*u + _mean + 0.445)); 282 } 283 284 if(v >= _u._ptrd.v_r) { 285 u = uniform_01<RealType>()(urng) - 0.5; 286 } else { 287 u = v/_u._ptrd.v_r - 0.93; 288 u = ((u < 0)? -0.5 : 0.5) - u; 289 v = uniform_01<RealType>()(urng) * _u._ptrd.v_r; 290 } 291 292 RealType us = 0.5 - abs(u); 293 if(us < 0.013 && v > us) { 294 continue; 295 } 296 297 RealType k = floor((2*_u._ptrd.a/us + _u._ptrd.b)*u+_mean+0.445); 298 v = v*_u._ptrd.inv_alpha/(_u._ptrd.a/(us*us) + _u._ptrd.b); 299 300 RealType log_sqrt_2pi = 0.91893853320467267; 301 302 if(k >= 10) { 303 if(log(v*_u._ptrd.smu) <= (k + 0.5)*log(_mean/k) 304 - _mean 305 - log_sqrt_2pi 306 + k 307 - (1/12. - (1/360. - 1/(1260.*k*k))/(k*k))/k) { 308 return static_cast<IntType>(k); 309 } 310 } else if(k >= 0) { 311 if(log(v) <= k*log(_mean) 312 - _mean 313 - log_factorial(static_cast<IntType>(k))) { 314 return static_cast<IntType>(k); 315 } 316 } 317 } 318 } 319 320 template<class URNG> invert(URNG & urng) const321 IntType invert(URNG& urng) const 322 { 323 RealType p = _u._exp_mean; 324 IntType x = 0; 325 RealType u = uniform_01<RealType>()(urng); 326 while(u > p) { 327 u = u - p; 328 ++x; 329 p = _mean * p / x; 330 } 331 return x; 332 } 333 334 RealType _mean; 335 336 union { 337 // for ptrd 338 struct { 339 RealType v_r; 340 RealType a; 341 RealType b; 342 RealType smu; 343 RealType inv_alpha; 344 } _ptrd; 345 // for inversion 346 RealType _exp_mean; 347 } _u; 348 349 /// @endcond 350 }; 351 352 } // namespace random 353 354 using random::poisson_distribution; 355 356 } // namespace boost 357 358 #include <boost/random/detail/enable_warnings.hpp> 359 360 #endif // BOOST_RANDOM_POISSON_DISTRIBUTION_HPP 361