]>
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, | |
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 | ||
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 |