1 // boost\math\distributions\poisson.hpp 2 3 // Copyright John Maddock 2006. 4 // Copyright Paul A. Bristow 2007. 5 6 // Use, modification and distribution are subject to the 7 // Boost Software License, Version 1.0. 8 // (See accompanying file LICENSE_1_0.txt 9 // or copy at http://www.boost.org/LICENSE_1_0.txt) 10 11 // Poisson distribution is a discrete probability distribution. 12 // It expresses the probability of a number (k) of 13 // events, occurrences, failures or arrivals occurring in a fixed time, 14 // assuming these events occur with a known average or mean rate (lambda) 15 // and are independent of the time since the last event. 16 // The distribution was discovered by Simeon-Denis Poisson (1781-1840). 17 18 // Parameter lambda is the mean number of events in the given time interval. 19 // The random variate k is the number of events, occurrences or arrivals. 20 // k argument may be integral, signed, or unsigned, or floating point. 21 // If necessary, it has already been promoted from an integral type. 22 23 // Note that the Poisson distribution 24 // (like others including the binomial, negative binomial & Bernoulli) 25 // is strictly defined as a discrete function: 26 // only integral values of k are envisaged. 27 // However because the method of calculation uses a continuous gamma function, 28 // it is convenient to treat it as if a continuous function, 29 // and permit non-integral values of k. 30 // To enforce the strict mathematical model, users should use floor or ceil functions 31 // on k outside this function to ensure that k is integral. 32 33 // See http://en.wikipedia.org/wiki/Poisson_distribution 34 // http://documents.wolfram.com/v5/Add-onsLinks/StandardPackages/Statistics/DiscreteDistributions.html 35 36 #ifndef BOOST_MATH_SPECIAL_POISSON_HPP 37 #define BOOST_MATH_SPECIAL_POISSON_HPP 38 39 #include <boost/math/distributions/fwd.hpp> 40 #include <boost/math/special_functions/gamma.hpp> // for incomplete gamma. gamma_q 41 #include <boost/math/special_functions/trunc.hpp> // for incomplete gamma. gamma_q 42 #include <boost/math/distributions/complement.hpp> // complements 43 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks 44 #include <boost/math/special_functions/fpclassify.hpp> // isnan. 45 #include <boost/math/special_functions/factorials.hpp> // factorials. 46 #include <boost/math/tools/roots.hpp> // for root finding. 47 #include <boost/math/distributions/detail/inv_discrete_quantile.hpp> 48 49 #include <utility> 50 51 namespace boost 52 { 53 namespace math 54 { 55 namespace poisson_detail 56 { 57 // Common error checking routines for Poisson distribution functions. 58 // These are convoluted, & apparently redundant, to try to ensure that 59 // checks are always performed, even if exceptions are not enabled. 60 61 template <class RealType, class Policy> check_mean(const char * function,const RealType & mean,RealType * result,const Policy & pol)62 inline bool check_mean(const char* function, const RealType& mean, RealType* result, const Policy& pol) 63 { 64 if(!(boost::math::isfinite)(mean) || (mean < 0)) 65 { 66 *result = policies::raise_domain_error<RealType>( 67 function, 68 "Mean argument is %1%, but must be >= 0 !", mean, pol); 69 return false; 70 } 71 return true; 72 } // bool check_mean 73 74 template <class RealType, class Policy> check_mean_NZ(const char * function,const RealType & mean,RealType * result,const Policy & pol)75 inline bool check_mean_NZ(const char* function, const RealType& mean, RealType* result, const Policy& pol) 76 { // mean == 0 is considered an error. 77 if( !(boost::math::isfinite)(mean) || (mean <= 0)) 78 { 79 *result = policies::raise_domain_error<RealType>( 80 function, 81 "Mean argument is %1%, but must be > 0 !", mean, pol); 82 return false; 83 } 84 return true; 85 } // bool check_mean_NZ 86 87 template <class RealType, class Policy> check_dist(const char * function,const RealType & mean,RealType * result,const Policy & pol)88 inline bool check_dist(const char* function, const RealType& mean, RealType* result, const Policy& pol) 89 { // Only one check, so this is redundant really but should be optimized away. 90 return check_mean_NZ(function, mean, result, pol); 91 } // bool check_dist 92 93 template <class RealType, class Policy> check_k(const char * function,const RealType & k,RealType * result,const Policy & pol)94 inline bool check_k(const char* function, const RealType& k, RealType* result, const Policy& pol) 95 { 96 if((k < 0) || !(boost::math::isfinite)(k)) 97 { 98 *result = policies::raise_domain_error<RealType>( 99 function, 100 "Number of events k argument is %1%, but must be >= 0 !", k, pol); 101 return false; 102 } 103 return true; 104 } // bool check_k 105 106 template <class RealType, class Policy> check_dist_and_k(const char * function,RealType mean,RealType k,RealType * result,const Policy & pol)107 inline bool check_dist_and_k(const char* function, RealType mean, RealType k, RealType* result, const Policy& pol) 108 { 109 if((check_dist(function, mean, result, pol) == false) || 110 (check_k(function, k, result, pol) == false)) 111 { 112 return false; 113 } 114 return true; 115 } // bool check_dist_and_k 116 117 template <class RealType, class Policy> check_prob(const char * function,const RealType & p,RealType * result,const Policy & pol)118 inline bool check_prob(const char* function, const RealType& p, RealType* result, const Policy& pol) 119 { // Check 0 <= p <= 1 120 if(!(boost::math::isfinite)(p) || (p < 0) || (p > 1)) 121 { 122 *result = policies::raise_domain_error<RealType>( 123 function, 124 "Probability argument is %1%, but must be >= 0 and <= 1 !", p, pol); 125 return false; 126 } 127 return true; 128 } // bool check_prob 129 130 template <class RealType, class Policy> check_dist_and_prob(const char * function,RealType mean,RealType p,RealType * result,const Policy & pol)131 inline bool check_dist_and_prob(const char* function, RealType mean, RealType p, RealType* result, const Policy& pol) 132 { 133 if((check_dist(function, mean, result, pol) == false) || 134 (check_prob(function, p, result, pol) == false)) 135 { 136 return false; 137 } 138 return true; 139 } // bool check_dist_and_prob 140 141 } // namespace poisson_detail 142 143 template <class RealType = double, class Policy = policies::policy<> > 144 class poisson_distribution 145 { 146 public: 147 typedef RealType value_type; 148 typedef Policy policy_type; 149 poisson_distribution(RealType l_mean=1)150 poisson_distribution(RealType l_mean = 1) : m_l(l_mean) // mean (lambda). 151 { // Expected mean number of events that occur during the given interval. 152 RealType r; 153 poisson_detail::check_dist( 154 "boost::math::poisson_distribution<%1%>::poisson_distribution", 155 m_l, 156 &r, Policy()); 157 } // poisson_distribution constructor. 158 mean() const159 RealType mean() const 160 { // Private data getter function. 161 return m_l; 162 } 163 private: 164 // Data member, initialized by constructor. 165 RealType m_l; // mean number of occurrences. 166 }; // template <class RealType, class Policy> class poisson_distribution 167 168 typedef poisson_distribution<double> poisson; // Reserved name of type double. 169 170 // Non-member functions to give properties of the distribution. 171 172 template <class RealType, class Policy> range(const poisson_distribution<RealType,Policy> &)173 inline const std::pair<RealType, RealType> range(const poisson_distribution<RealType, Policy>& /* dist */) 174 { // Range of permissible values for random variable k. 175 using boost::math::tools::max_value; 176 return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // Max integer? 177 } 178 179 template <class RealType, class Policy> support(const poisson_distribution<RealType,Policy> &)180 inline const std::pair<RealType, RealType> support(const poisson_distribution<RealType, Policy>& /* dist */) 181 { // Range of supported values for random variable k. 182 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. 183 using boost::math::tools::max_value; 184 return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); 185 } 186 187 template <class RealType, class Policy> mean(const poisson_distribution<RealType,Policy> & dist)188 inline RealType mean(const poisson_distribution<RealType, Policy>& dist) 189 { // Mean of poisson distribution = lambda. 190 return dist.mean(); 191 } // mean 192 193 template <class RealType, class Policy> mode(const poisson_distribution<RealType,Policy> & dist)194 inline RealType mode(const poisson_distribution<RealType, Policy>& dist) 195 { // mode. 196 BOOST_MATH_STD_USING // ADL of std functions. 197 return floor(dist.mean()); 198 } 199 200 //template <class RealType, class Policy> 201 //inline RealType median(const poisson_distribution<RealType, Policy>& dist) 202 //{ // median = approximately lambda + 1/3 - 0.2/lambda 203 // RealType l = dist.mean(); 204 // return dist.mean() + static_cast<RealType>(0.3333333333333333333333333333333333333333333333) 205 // - static_cast<RealType>(0.2) / l; 206 //} // BUT this formula appears to be out-by-one compared to quantile(half) 207 // Query posted on Wikipedia. 208 // Now implemented via quantile(half) in derived accessors. 209 210 template <class RealType, class Policy> variance(const poisson_distribution<RealType,Policy> & dist)211 inline RealType variance(const poisson_distribution<RealType, Policy>& dist) 212 { // variance. 213 return dist.mean(); 214 } 215 216 // RealType standard_deviation(const poisson_distribution<RealType, Policy>& dist) 217 // standard_deviation provided by derived accessors. 218 219 template <class RealType, class Policy> skewness(const poisson_distribution<RealType,Policy> & dist)220 inline RealType skewness(const poisson_distribution<RealType, Policy>& dist) 221 { // skewness = sqrt(l). 222 BOOST_MATH_STD_USING // ADL of std functions. 223 return 1 / sqrt(dist.mean()); 224 } 225 226 template <class RealType, class Policy> kurtosis_excess(const poisson_distribution<RealType,Policy> & dist)227 inline RealType kurtosis_excess(const poisson_distribution<RealType, Policy>& dist) 228 { // skewness = sqrt(l). 229 return 1 / dist.mean(); // kurtosis_excess 1/mean from Wiki & MathWorld eq 31. 230 // http://mathworld.wolfram.com/Kurtosis.html explains that the kurtosis excess 231 // is more convenient because the kurtosis excess of a normal distribution is zero 232 // whereas the true kurtosis is 3. 233 } // RealType kurtosis_excess 234 235 template <class RealType, class Policy> kurtosis(const poisson_distribution<RealType,Policy> & dist)236 inline RealType kurtosis(const poisson_distribution<RealType, Policy>& dist) 237 { // kurtosis is 4th moment about the mean = u4 / sd ^ 4 238 // http://en.wikipedia.org/wiki/Kurtosis 239 // kurtosis can range from -2 (flat top) to +infinity (sharp peak & heavy tails). 240 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm 241 return 3 + 1 / dist.mean(); // NIST. 242 // http://mathworld.wolfram.com/Kurtosis.html explains that the kurtosis excess 243 // is more convenient because the kurtosis excess of a normal distribution is zero 244 // whereas the true kurtosis is 3. 245 } // RealType kurtosis 246 247 template <class RealType, class Policy> 248 RealType pdf(const poisson_distribution<RealType, Policy>& dist, const RealType& k) 249 { // Probability Density/Mass Function. 250 // Probability that there are EXACTLY k occurrences (or arrivals). 251 BOOST_FPU_EXCEPTION_GUARD 252 253 BOOST_MATH_STD_USING // for ADL of std functions. 254 255 RealType mean = dist.mean(); 256 // Error check: 257 RealType result = 0; 258 if(false == poisson_detail::check_dist_and_k( 259 "boost::math::pdf(const poisson_distribution<%1%>&, %1%)", 260 mean, 261 k, 262 &result, Policy())) 263 { 264 return result; 265 } 266 267 // Special case of mean zero, regardless of the number of events k. 268 if (mean == 0) 269 { // Probability for any k is zero. 270 return 0; 271 } 272 if (k == 0) 273 { // mean ^ k = 1, and k! = 1, so can simplify. 274 return exp(-mean); 275 } 276 return boost::math::gamma_p_derivative(k+1, mean, Policy()); 277 } // pdf 278 279 template <class RealType, class Policy> 280 RealType cdf(const poisson_distribution<RealType, Policy>& dist, const RealType& k) 281 { // Cumulative Distribution Function Poisson. 282 // The random variate k is the number of occurrences(or arrivals) 283 // k argument may be integral, signed, or unsigned, or floating point. 284 // If necessary, it has already been promoted from an integral type. 285 // Returns the sum of the terms 0 through k of the Poisson Probability Density or Mass (pdf). 286 287 // But note that the Poisson distribution 288 // (like others including the binomial, negative binomial & Bernoulli) 289 // is strictly defined as a discrete function: only integral values of k are envisaged. 290 // However because of the method of calculation using a continuous gamma function, 291 // it is convenient to treat it as if it is a continuous function 292 // and permit non-integral values of k. 293 // To enforce the strict mathematical model, users should use floor or ceil functions 294 // outside this function to ensure that k is integral. 295 296 // The terms are not summed directly (at least for larger k) 297 // instead the incomplete gamma integral is employed, 298 299 BOOST_MATH_STD_USING // for ADL of std function exp. 300 301 RealType mean = dist.mean(); 302 // Error checks: 303 RealType result = 0; 304 if(false == poisson_detail::check_dist_and_k( 305 "boost::math::cdf(const poisson_distribution<%1%>&, %1%)", 306 mean, 307 k, 308 &result, Policy())) 309 { 310 return result; 311 } 312 // Special cases: 313 if (mean == 0) 314 { // Probability for any k is zero. 315 return 0; 316 } 317 if (k == 0) 318 { // return pdf(dist, static_cast<RealType>(0)); 319 // but mean (and k) have already been checked, 320 // so this avoids unnecessary repeated checks. 321 return exp(-mean); 322 } 323 // For small integral k could use a finite sum - 324 // it's cheaper than the gamma function. 325 // BUT this is now done efficiently by gamma_q function. 326 // Calculate poisson cdf using the gamma_q function. 327 return gamma_q(k+1, mean, Policy()); 328 } // binomial cdf 329 330 template <class RealType, class Policy> 331 RealType cdf(const complemented2_type<poisson_distribution<RealType, Policy>, RealType>& c) 332 { // Complemented Cumulative Distribution Function Poisson 333 // The random variate k is the number of events, occurrences or arrivals. 334 // k argument may be integral, signed, or unsigned, or floating point. 335 // If necessary, it has already been promoted from an integral type. 336 // But note that the Poisson distribution 337 // (like others including the binomial, negative binomial & Bernoulli) 338 // is strictly defined as a discrete function: only integral values of k are envisaged. 339 // However because of the method of calculation using a continuous gamma function, 340 // it is convenient to treat it as is it is a continuous function 341 // and permit non-integral values of k. 342 // To enforce the strict mathematical model, users should use floor or ceil functions 343 // outside this function to ensure that k is integral. 344 345 // Returns the sum of the terms k+1 through inf of the Poisson Probability Density/Mass (pdf). 346 // The terms are not summed directly (at least for larger k) 347 // instead the incomplete gamma integral is employed, 348 349 RealType const& k = c.param; 350 poisson_distribution<RealType, Policy> const& dist = c.dist; 351 352 RealType mean = dist.mean(); 353 354 // Error checks: 355 RealType result = 0; 356 if(false == poisson_detail::check_dist_and_k( 357 "boost::math::cdf(const poisson_distribution<%1%>&, %1%)", 358 mean, 359 k, 360 &result, Policy())) 361 { 362 return result; 363 } 364 // Special case of mean, regardless of the number of events k. 365 if (mean == 0) 366 { // Probability for any k is unity, complement of zero. 367 return 1; 368 } 369 if (k == 0) 370 { // Avoid repeated checks on k and mean in gamma_p. 371 return -boost::math::expm1(-mean, Policy()); 372 } 373 // Unlike un-complemented cdf (sum from 0 to k), 374 // can't use finite sum from k+1 to infinity for small integral k, 375 // anyway it is now done efficiently by gamma_p. 376 return gamma_p(k + 1, mean, Policy()); // Calculate Poisson cdf using the gamma_p function. 377 // CCDF = gamma_p(k+1, lambda) 378 } // poisson ccdf 379 380 template <class RealType, class Policy> quantile(const poisson_distribution<RealType,Policy> & dist,const RealType & p)381 inline RealType quantile(const poisson_distribution<RealType, Policy>& dist, const RealType& p) 382 { // Quantile (or Percent Point) Poisson function. 383 // Return the number of expected events k for a given probability p. 384 static const char* function = "boost::math::quantile(const poisson_distribution<%1%>&, %1%)"; 385 RealType result = 0; // of Argument checks: 386 if(false == poisson_detail::check_prob( 387 function, 388 p, 389 &result, Policy())) 390 { 391 return result; 392 } 393 // Special case: 394 if (dist.mean() == 0) 395 { // if mean = 0 then p = 0, so k can be anything? 396 if (false == poisson_detail::check_mean_NZ( 397 function, 398 dist.mean(), 399 &result, Policy())) 400 { 401 return result; 402 } 403 } 404 if(p == 0) 405 { 406 return 0; // Exact result regardless of discrete-quantile Policy 407 } 408 if(p == 1) 409 { 410 return policies::raise_overflow_error<RealType>(function, 0, Policy()); 411 } 412 typedef typename Policy::discrete_quantile_type discrete_type; 413 boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); 414 RealType guess, factor = 8; 415 RealType z = dist.mean(); 416 if(z < 1) 417 guess = z; 418 else 419 guess = boost::math::detail::inverse_poisson_cornish_fisher(z, p, RealType(1-p), Policy()); 420 if(z > 5) 421 { 422 if(z > 1000) 423 factor = 1.01f; 424 else if(z > 50) 425 factor = 1.1f; 426 else if(guess > 10) 427 factor = 1.25f; 428 else 429 factor = 2; 430 if(guess < 1.1) 431 factor = 8; 432 } 433 434 return detail::inverse_discrete_quantile( 435 dist, 436 p, 437 false, 438 guess, 439 factor, 440 RealType(1), 441 discrete_type(), 442 max_iter); 443 } // quantile 444 445 template <class RealType, class Policy> quantile(const complemented2_type<poisson_distribution<RealType,Policy>,RealType> & c)446 inline RealType quantile(const complemented2_type<poisson_distribution<RealType, Policy>, RealType>& c) 447 { // Quantile (or Percent Point) of Poisson function. 448 // Return the number of expected events k for a given 449 // complement of the probability q. 450 // 451 // Error checks: 452 static const char* function = "boost::math::quantile(complement(const poisson_distribution<%1%>&, %1%))"; 453 RealType q = c.param; 454 const poisson_distribution<RealType, Policy>& dist = c.dist; 455 RealType result = 0; // of argument checks. 456 if(false == poisson_detail::check_prob( 457 function, 458 q, 459 &result, Policy())) 460 { 461 return result; 462 } 463 // Special case: 464 if (dist.mean() == 0) 465 { // if mean = 0 then p = 0, so k can be anything? 466 if (false == poisson_detail::check_mean_NZ( 467 function, 468 dist.mean(), 469 &result, Policy())) 470 { 471 return result; 472 } 473 } 474 if(q == 0) 475 { 476 return policies::raise_overflow_error<RealType>(function, 0, Policy()); 477 } 478 if(q == 1) 479 { 480 return 0; // Exact result regardless of discrete-quantile Policy 481 } 482 typedef typename Policy::discrete_quantile_type discrete_type; 483 boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); 484 RealType guess, factor = 8; 485 RealType z = dist.mean(); 486 if(z < 1) 487 guess = z; 488 else 489 guess = boost::math::detail::inverse_poisson_cornish_fisher(z, RealType(1-q), q, Policy()); 490 if(z > 5) 491 { 492 if(z > 1000) 493 factor = 1.01f; 494 else if(z > 50) 495 factor = 1.1f; 496 else if(guess > 10) 497 factor = 1.25f; 498 else 499 factor = 2; 500 if(guess < 1.1) 501 factor = 8; 502 } 503 504 return detail::inverse_discrete_quantile( 505 dist, 506 q, 507 true, 508 guess, 509 factor, 510 RealType(1), 511 discrete_type(), 512 max_iter); 513 } // quantile complement. 514 515 } // namespace math 516 } // namespace boost 517 518 // This include must be at the end, *after* the accessors 519 // for this distribution have been defined, in order to 520 // keep compilers that support two-phase lookup happy. 521 #include <boost/math/distributions/detail/derived_accessors.hpp> 522 #include <boost/math/distributions/detail/inv_discrete_quantile.hpp> 523 524 #endif // BOOST_MATH_SPECIAL_POISSON_HPP 525 526 527 528