ceph/src/boost/libs/math/include/boost/math/distributions/poisson.hpp

   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 continous 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>
  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>
  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>
  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>
  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>
 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>
 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>
 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
 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
 159       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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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/Curtosis
 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 continous 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 continous 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>
 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>
 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