]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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, | |
f67539c2 | 28 | // it is convenient to treat it as if a continuous function, |
7c673cae FG |
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 | ||
1e59de90 TL |
170 | #ifdef __cpp_deduction_guides |
171 | template <class RealType> | |
172 | poisson_distribution(RealType)->poisson_distribution<typename boost::math::tools::promote_args<RealType>::type>; | |
173 | #endif | |
174 | ||
7c673cae FG |
175 | // Non-member functions to give properties of the distribution. |
176 | ||
177 | template <class RealType, class Policy> | |
178 | inline const std::pair<RealType, RealType> range(const poisson_distribution<RealType, Policy>& /* dist */) | |
179 | { // Range of permissible values for random variable k. | |
180 | using boost::math::tools::max_value; | |
181 | return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // Max integer? | |
182 | } | |
183 | ||
184 | template <class RealType, class Policy> | |
185 | inline const std::pair<RealType, RealType> support(const poisson_distribution<RealType, Policy>& /* dist */) | |
186 | { // Range of supported values for random variable k. | |
187 | // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. | |
188 | using boost::math::tools::max_value; | |
189 | return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); | |
190 | } | |
191 | ||
192 | template <class RealType, class Policy> | |
193 | inline RealType mean(const poisson_distribution<RealType, Policy>& dist) | |
194 | { // Mean of poisson distribution = lambda. | |
195 | return dist.mean(); | |
196 | } // mean | |
197 | ||
198 | template <class RealType, class Policy> | |
199 | inline RealType mode(const poisson_distribution<RealType, Policy>& dist) | |
200 | { // mode. | |
201 | BOOST_MATH_STD_USING // ADL of std functions. | |
202 | return floor(dist.mean()); | |
203 | } | |
204 | ||
205 | //template <class RealType, class Policy> | |
206 | //inline RealType median(const poisson_distribution<RealType, Policy>& dist) | |
207 | //{ // median = approximately lambda + 1/3 - 0.2/lambda | |
208 | // RealType l = dist.mean(); | |
209 | // return dist.mean() + static_cast<RealType>(0.3333333333333333333333333333333333333333333333) | |
210 | // - static_cast<RealType>(0.2) / l; | |
211 | //} // BUT this formula appears to be out-by-one compared to quantile(half) | |
212 | // Query posted on Wikipedia. | |
213 | // Now implemented via quantile(half) in derived accessors. | |
214 | ||
215 | template <class RealType, class Policy> | |
216 | inline RealType variance(const poisson_distribution<RealType, Policy>& dist) | |
217 | { // variance. | |
218 | return dist.mean(); | |
219 | } | |
220 | ||
221 | // RealType standard_deviation(const poisson_distribution<RealType, Policy>& dist) | |
222 | // standard_deviation provided by derived accessors. | |
223 | ||
224 | template <class RealType, class Policy> | |
225 | inline RealType skewness(const poisson_distribution<RealType, Policy>& dist) | |
226 | { // skewness = sqrt(l). | |
227 | BOOST_MATH_STD_USING // ADL of std functions. | |
228 | return 1 / sqrt(dist.mean()); | |
229 | } | |
230 | ||
231 | template <class RealType, class Policy> | |
232 | inline RealType kurtosis_excess(const poisson_distribution<RealType, Policy>& dist) | |
233 | { // skewness = sqrt(l). | |
234 | return 1 / dist.mean(); // kurtosis_excess 1/mean from Wiki & MathWorld eq 31. | |
235 | // http://mathworld.wolfram.com/Kurtosis.html explains that the kurtosis excess | |
236 | // is more convenient because the kurtosis excess of a normal distribution is zero | |
237 | // whereas the true kurtosis is 3. | |
238 | } // RealType kurtosis_excess | |
239 | ||
240 | template <class RealType, class Policy> | |
241 | inline RealType kurtosis(const poisson_distribution<RealType, Policy>& dist) | |
242 | { // kurtosis is 4th moment about the mean = u4 / sd ^ 4 | |
f67539c2 | 243 | // http://en.wikipedia.org/wiki/Kurtosis |
7c673cae FG |
244 | // kurtosis can range from -2 (flat top) to +infinity (sharp peak & heavy tails). |
245 | // http://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm | |
246 | return 3 + 1 / dist.mean(); // NIST. | |
247 | // http://mathworld.wolfram.com/Kurtosis.html explains that the kurtosis excess | |
248 | // is more convenient because the kurtosis excess of a normal distribution is zero | |
249 | // whereas the true kurtosis is 3. | |
250 | } // RealType kurtosis | |
251 | ||
252 | template <class RealType, class Policy> | |
253 | RealType pdf(const poisson_distribution<RealType, Policy>& dist, const RealType& k) | |
254 | { // Probability Density/Mass Function. | |
255 | // Probability that there are EXACTLY k occurrences (or arrivals). | |
256 | BOOST_FPU_EXCEPTION_GUARD | |
257 | ||
258 | BOOST_MATH_STD_USING // for ADL of std functions. | |
259 | ||
260 | RealType mean = dist.mean(); | |
261 | // Error check: | |
262 | RealType result = 0; | |
263 | if(false == poisson_detail::check_dist_and_k( | |
264 | "boost::math::pdf(const poisson_distribution<%1%>&, %1%)", | |
265 | mean, | |
266 | k, | |
267 | &result, Policy())) | |
268 | { | |
269 | return result; | |
270 | } | |
271 | ||
272 | // Special case of mean zero, regardless of the number of events k. | |
273 | if (mean == 0) | |
274 | { // Probability for any k is zero. | |
275 | return 0; | |
276 | } | |
277 | if (k == 0) | |
278 | { // mean ^ k = 1, and k! = 1, so can simplify. | |
279 | return exp(-mean); | |
280 | } | |
281 | return boost::math::gamma_p_derivative(k+1, mean, Policy()); | |
282 | ||
283 | ||
284 | template <class RealType, class Policy> | |
285 | RealType cdf(const poisson_distribution<RealType, Policy>& dist, const RealType& k) | |
286 | { // Cumulative Distribution Function Poisson. | |
287 | // The random variate k is the number of occurrences(or arrivals) | |
288 | // k argument may be integral, signed, or unsigned, or floating point. | |
289 | // If necessary, it has already been promoted from an integral type. | |
290 | // Returns the sum of the terms 0 through k of the Poisson Probability Density or Mass (pdf). | |
291 | ||
292 | // But note that the Poisson distribution | |
293 | // (like others including the binomial, negative binomial & Bernoulli) | |
294 | // is strictly defined as a discrete function: only integral values of k are envisaged. | |
295 | // However because of the method of calculation using a continuous gamma function, | |
f67539c2 | 296 | // it is convenient to treat it as if it is a continuous function |
7c673cae FG |
297 | // and permit non-integral values of k. |
298 | // To enforce the strict mathematical model, users should use floor or ceil functions | |
299 | // outside this function to ensure that k is integral. | |
300 | ||
301 | // The terms are not summed directly (at least for larger k) | |
302 | // instead the incomplete gamma integral is employed, | |
303 | ||
304 | BOOST_MATH_STD_USING // for ADL of std function exp. | |
305 | ||
306 | RealType mean = dist.mean(); | |
307 | // Error checks: | |
308 | RealType result = 0; | |
309 | if(false == poisson_detail::check_dist_and_k( | |
310 | "boost::math::cdf(const poisson_distribution<%1%>&, %1%)", | |
311 | mean, | |
312 | k, | |
313 | &result, Policy())) | |
314 | { | |
315 | return result; | |
316 | } | |
317 | // Special cases: | |
318 | if (mean == 0) | |
319 | { // Probability for any k is zero. | |
320 | return 0; | |
321 | } | |
322 | if (k == 0) | |
323 | { // return pdf(dist, static_cast<RealType>(0)); | |
324 | // but mean (and k) have already been checked, | |
325 | // so this avoids unnecessary repeated checks. | |
326 | return exp(-mean); | |
327 | } | |
328 | // For small integral k could use a finite sum - | |
329 | // it's cheaper than the gamma function. | |
330 | // BUT this is now done efficiently by gamma_q function. | |
331 | // Calculate poisson cdf using the gamma_q function. | |
332 | return gamma_q(k+1, mean, Policy()); | |
333 | } // binomial cdf | |
334 | ||
335 | template <class RealType, class Policy> | |
336 | RealType cdf(const complemented2_type<poisson_distribution<RealType, Policy>, RealType>& c) | |
337 | { // Complemented Cumulative Distribution Function Poisson | |
338 | // The random variate k is the number of events, occurrences or arrivals. | |
339 | // k argument may be integral, signed, or unsigned, or floating point. | |
340 | // If necessary, it has already been promoted from an integral type. | |
341 | // But note that the Poisson distribution | |
342 | // (like others including the binomial, negative binomial & Bernoulli) | |
343 | // is strictly defined as a discrete function: only integral values of k are envisaged. | |
344 | // However because of the method of calculation using a continuous gamma function, | |
f67539c2 | 345 | // it is convenient to treat it as is it is a continuous function |
7c673cae FG |
346 | // and permit non-integral values of k. |
347 | // To enforce the strict mathematical model, users should use floor or ceil functions | |
348 | // outside this function to ensure that k is integral. | |
349 | ||
350 | // Returns the sum of the terms k+1 through inf of the Poisson Probability Density/Mass (pdf). | |
351 | // The terms are not summed directly (at least for larger k) | |
352 | // instead the incomplete gamma integral is employed, | |
353 | ||
354 | RealType const& k = c.param; | |
355 | poisson_distribution<RealType, Policy> const& dist = c.dist; | |
356 | ||
357 | RealType mean = dist.mean(); | |
358 | ||
359 | // Error checks: | |
360 | RealType result = 0; | |
361 | if(false == poisson_detail::check_dist_and_k( | |
362 | "boost::math::cdf(const poisson_distribution<%1%>&, %1%)", | |
363 | mean, | |
364 | k, | |
365 | &result, Policy())) | |
366 | { | |
367 | return result; | |
368 | } | |
369 | // Special case of mean, regardless of the number of events k. | |
370 | if (mean == 0) | |
371 | { // Probability for any k is unity, complement of zero. | |
372 | return 1; | |
373 | } | |
374 | if (k == 0) | |
375 | { // Avoid repeated checks on k and mean in gamma_p. | |
376 | return -boost::math::expm1(-mean, Policy()); | |
377 | } | |
378 | // Unlike un-complemented cdf (sum from 0 to k), | |
379 | // can't use finite sum from k+1 to infinity for small integral k, | |
380 | // anyway it is now done efficiently by gamma_p. | |
381 | return gamma_p(k + 1, mean, Policy()); // Calculate Poisson cdf using the gamma_p function. | |
382 | // CCDF = gamma_p(k+1, lambda) | |
383 | } // poisson ccdf | |
384 | ||
385 | template <class RealType, class Policy> | |
386 | inline RealType quantile(const poisson_distribution<RealType, Policy>& dist, const RealType& p) | |
387 | { // Quantile (or Percent Point) Poisson function. | |
388 | // Return the number of expected events k for a given probability p. | |
389 | static const char* function = "boost::math::quantile(const poisson_distribution<%1%>&, %1%)"; | |
390 | RealType result = 0; // of Argument checks: | |
391 | if(false == poisson_detail::check_prob( | |
392 | function, | |
393 | p, | |
394 | &result, Policy())) | |
395 | { | |
396 | return result; | |
397 | } | |
398 | // Special case: | |
399 | if (dist.mean() == 0) | |
400 | { // if mean = 0 then p = 0, so k can be anything? | |
401 | if (false == poisson_detail::check_mean_NZ( | |
402 | function, | |
403 | dist.mean(), | |
404 | &result, Policy())) | |
405 | { | |
406 | return result; | |
407 | } | |
408 | } | |
409 | if(p == 0) | |
410 | { | |
411 | return 0; // Exact result regardless of discrete-quantile Policy | |
412 | } | |
413 | if(p == 1) | |
414 | { | |
415 | return policies::raise_overflow_error<RealType>(function, 0, Policy()); | |
416 | } | |
417 | typedef typename Policy::discrete_quantile_type discrete_type; | |
1e59de90 | 418 | std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); |
7c673cae FG |
419 | RealType guess, factor = 8; |
420 | RealType z = dist.mean(); | |
421 | if(z < 1) | |
422 | guess = z; | |
423 | else | |
424 | guess = boost::math::detail::inverse_poisson_cornish_fisher(z, p, RealType(1-p), Policy()); | |
425 | if(z > 5) | |
426 | { | |
427 | if(z > 1000) | |
428 | factor = 1.01f; | |
429 | else if(z > 50) | |
430 | factor = 1.1f; | |
431 | else if(guess > 10) | |
432 | factor = 1.25f; | |
433 | else | |
434 | factor = 2; | |
435 | if(guess < 1.1) | |
436 | factor = 8; | |
437 | } | |
438 | ||
439 | return detail::inverse_discrete_quantile( | |
440 | dist, | |
441 | p, | |
442 | false, | |
443 | guess, | |
444 | factor, | |
445 | RealType(1), | |
446 | discrete_type(), | |
447 | max_iter); | |
448 | } // quantile | |
449 | ||
450 | template <class RealType, class Policy> | |
451 | inline RealType quantile(const complemented2_type<poisson_distribution<RealType, Policy>, RealType>& c) | |
452 | { // Quantile (or Percent Point) of Poisson function. | |
453 | // Return the number of expected events k for a given | |
454 | // complement of the probability q. | |
455 | // | |
456 | // Error checks: | |
457 | static const char* function = "boost::math::quantile(complement(const poisson_distribution<%1%>&, %1%))"; | |
458 | RealType q = c.param; | |
459 | const poisson_distribution<RealType, Policy>& dist = c.dist; | |
460 | RealType result = 0; // of argument checks. | |
461 | if(false == poisson_detail::check_prob( | |
462 | function, | |
463 | q, | |
464 | &result, Policy())) | |
465 | { | |
466 | return result; | |
467 | } | |
468 | // Special case: | |
469 | if (dist.mean() == 0) | |
470 | { // if mean = 0 then p = 0, so k can be anything? | |
471 | if (false == poisson_detail::check_mean_NZ( | |
472 | function, | |
473 | dist.mean(), | |
474 | &result, Policy())) | |
475 | { | |
476 | return result; | |
477 | } | |
478 | } | |
479 | if(q == 0) | |
480 | { | |
481 | return policies::raise_overflow_error<RealType>(function, 0, Policy()); | |
482 | } | |
483 | if(q == 1) | |
484 | { | |
485 | return 0; // Exact result regardless of discrete-quantile Policy | |
486 | } | |
487 | typedef typename Policy::discrete_quantile_type discrete_type; | |
1e59de90 | 488 | std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); |
7c673cae FG |
489 | RealType guess, factor = 8; |
490 | RealType z = dist.mean(); | |
491 | if(z < 1) | |
492 | guess = z; | |
493 | else | |
494 | guess = boost::math::detail::inverse_poisson_cornish_fisher(z, RealType(1-q), q, Policy()); | |
495 | if(z > 5) | |
496 | { | |
497 | if(z > 1000) | |
498 | factor = 1.01f; | |
499 | else if(z > 50) | |
500 | factor = 1.1f; | |
501 | else if(guess > 10) | |
502 | factor = 1.25f; | |
503 | else | |
504 | factor = 2; | |
505 | if(guess < 1.1) | |
506 | factor = 8; | |
507 | } | |
508 | ||
509 | return detail::inverse_discrete_quantile( | |
510 | dist, | |
511 | q, | |
512 | true, | |
513 | guess, | |
514 | factor, | |
515 | RealType(1), | |
516 | discrete_type(), | |
517 | max_iter); | |
518 | } // quantile complement. | |
519 | ||
520 | } // namespace math | |
521 | } // namespace boost | |
522 | ||
523 | // This include must be at the end, *after* the accessors | |
524 | // for this distribution have been defined, in order to | |
525 | // keep compilers that support two-phase lookup happy. | |
526 | #include <boost/math/distributions/detail/derived_accessors.hpp> | |
527 | #include <boost/math/distributions/detail/inv_discrete_quantile.hpp> | |
528 | ||
529 | #endif // BOOST_MATH_SPECIAL_POISSON_HPP | |
530 | ||
531 | ||
532 |