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

source code of boost/libs/math/include/boost/math/distributions/poisson.hpp