| // boost\math\distributions\poisson.hpp |
| |
| // Copyright John Maddock 2006. |
| // Copyright Paul A. Bristow 2007. |
| |
| // Use, modification and distribution are subject to the |
| // Boost Software License, Version 1.0. |
| // (See accompanying file LICENSE_1_0.txt |
| // or copy at http://www.boost.org/LICENSE_1_0.txt) |
| |
| // Poisson distribution is a discrete probability distribution. |
| // It expresses the probability of a number (k) of |
| // events, occurrences, failures or arrivals occurring in a fixed time, |
| // assuming these events occur with a known average or mean rate (lambda) |
| // and are independent of the time since the last event. |
| // The distribution was discovered by Simeon-Denis Poisson (1781-1840). |
| |
| // Parameter lambda is the mean number of events in the given time interval. |
| // The random variate k is the number of events, occurrences or arrivals. |
| // k argument may be integral, signed, or unsigned, or floating point. |
| // If necessary, it has already been promoted from an integral type. |
| |
| // Note that the Poisson distribution |
| // (like others including the binomial, negative binomial & Bernoulli) |
| // is strictly defined as a discrete function: |
| // only integral values of k are envisaged. |
| // However because the method of calculation uses a continuous gamma function, |
| // it is convenient to treat it as if a continous function, |
| // and permit non-integral values of k. |
| // To enforce the strict mathematical model, users should use floor or ceil functions |
| // on k outside this function to ensure that k is integral. |
| |
| // See http://en.wikipedia.org/wiki/Poisson_distribution |
| // http://documents.wolfram.com/v5/Add-onsLinks/StandardPackages/Statistics/DiscreteDistributions.html |
| |
| #ifndef BOOST_MATH_SPECIAL_POISSON_HPP |
| #define BOOST_MATH_SPECIAL_POISSON_HPP |
| |
| #include <boost/math/distributions/fwd.hpp> |
| #include <boost/math/special_functions/gamma.hpp> // for incomplete gamma. gamma_q |
| #include <boost/math/special_functions/trunc.hpp> // for incomplete gamma. gamma_q |
| #include <boost/math/distributions/complement.hpp> // complements |
| #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks |
| #include <boost/math/special_functions/fpclassify.hpp> // isnan. |
| #include <boost/math/special_functions/factorials.hpp> // factorials. |
| #include <boost/math/tools/roots.hpp> // for root finding. |
| #include <boost/math/distributions/detail/inv_discrete_quantile.hpp> |
| |
| #include <utility> |
| |
| namespace boost |
| { |
| namespace math |
| { |
| namespace detail{ |
| template <class Dist> |
| inline typename Dist::value_type |
| inverse_discrete_quantile( |
| const Dist& dist, |
| const typename Dist::value_type& p, |
| const typename Dist::value_type& guess, |
| const typename Dist::value_type& multiplier, |
| const typename Dist::value_type& adder, |
| const policies::discrete_quantile<policies::integer_round_nearest>&, |
| boost::uintmax_t& max_iter); |
| template <class Dist> |
| inline typename Dist::value_type |
| inverse_discrete_quantile( |
| const Dist& dist, |
| const typename Dist::value_type& p, |
| const typename Dist::value_type& guess, |
| const typename Dist::value_type& multiplier, |
| const typename Dist::value_type& adder, |
| const policies::discrete_quantile<policies::integer_round_up>&, |
| boost::uintmax_t& max_iter); |
| template <class Dist> |
| inline typename Dist::value_type |
| inverse_discrete_quantile( |
| const Dist& dist, |
| const typename Dist::value_type& p, |
| const typename Dist::value_type& guess, |
| const typename Dist::value_type& multiplier, |
| const typename Dist::value_type& adder, |
| const policies::discrete_quantile<policies::integer_round_down>&, |
| boost::uintmax_t& max_iter); |
| template <class Dist> |
| inline typename Dist::value_type |
| inverse_discrete_quantile( |
| const Dist& dist, |
| const typename Dist::value_type& p, |
| const typename Dist::value_type& guess, |
| const typename Dist::value_type& multiplier, |
| const typename Dist::value_type& adder, |
| const policies::discrete_quantile<policies::integer_round_outwards>&, |
| boost::uintmax_t& max_iter); |
| template <class Dist> |
| inline typename Dist::value_type |
| inverse_discrete_quantile( |
| const Dist& dist, |
| const typename Dist::value_type& p, |
| const typename Dist::value_type& guess, |
| const typename Dist::value_type& multiplier, |
| const typename Dist::value_type& adder, |
| const policies::discrete_quantile<policies::integer_round_inwards>&, |
| boost::uintmax_t& max_iter); |
| template <class Dist> |
| inline typename Dist::value_type |
| inverse_discrete_quantile( |
| const Dist& dist, |
| const typename Dist::value_type& p, |
| const typename Dist::value_type& guess, |
| const typename Dist::value_type& multiplier, |
| const typename Dist::value_type& adder, |
| const policies::discrete_quantile<policies::real>&, |
| boost::uintmax_t& max_iter); |
| } |
| namespace poisson_detail |
| { |
| // Common error checking routines for Poisson distribution functions. |
| // These are convoluted, & apparently redundant, to try to ensure that |
| // checks are always performed, even if exceptions are not enabled. |
| |
| template <class RealType, class Policy> |
| inline bool check_mean(const char* function, const RealType& mean, RealType* result, const Policy& pol) |
| { |
| if(!(boost::math::isfinite)(mean) || (mean < 0)) |
| { |
| *result = policies::raise_domain_error<RealType>( |
| function, |
| "Mean argument is %1%, but must be >= 0 !", mean, pol); |
| return false; |
| } |
| return true; |
| } // bool check_mean |
| |
| template <class RealType, class Policy> |
| inline bool check_mean_NZ(const char* function, const RealType& mean, RealType* result, const Policy& pol) |
| { // mean == 0 is considered an error. |
| if( !(boost::math::isfinite)(mean) || (mean <= 0)) |
| { |
| *result = policies::raise_domain_error<RealType>( |
| function, |
| "Mean argument is %1%, but must be > 0 !", mean, pol); |
| return false; |
| } |
| return true; |
| } // bool check_mean_NZ |
| |
| template <class RealType, class Policy> |
| inline bool check_dist(const char* function, const RealType& mean, RealType* result, const Policy& pol) |
| { // Only one check, so this is redundant really but should be optimized away. |
| return check_mean_NZ(function, mean, result, pol); |
| } // bool check_dist |
| |
| template <class RealType, class Policy> |
| inline bool check_k(const char* function, const RealType& k, RealType* result, const Policy& pol) |
| { |
| if((k < 0) || !(boost::math::isfinite)(k)) |
| { |
| *result = policies::raise_domain_error<RealType>( |
| function, |
| "Number of events k argument is %1%, but must be >= 0 !", k, pol); |
| return false; |
| } |
| return true; |
| } // bool check_k |
| |
| template <class RealType, class Policy> |
| inline bool check_dist_and_k(const char* function, RealType mean, RealType k, RealType* result, const Policy& pol) |
| { |
| if((check_dist(function, mean, result, pol) == false) || |
| (check_k(function, k, result, pol) == false)) |
| { |
| return false; |
| } |
| return true; |
| } // bool check_dist_and_k |
| |
| template <class RealType, class Policy> |
| inline bool check_prob(const char* function, const RealType& p, RealType* result, const Policy& pol) |
| { // Check 0 <= p <= 1 |
| if(!(boost::math::isfinite)(p) || (p < 0) || (p > 1)) |
| { |
| *result = policies::raise_domain_error<RealType>( |
| function, |
| "Probability argument is %1%, but must be >= 0 and <= 1 !", p, pol); |
| return false; |
| } |
| return true; |
| } // bool check_prob |
| |
| template <class RealType, class Policy> |
| inline bool check_dist_and_prob(const char* function, RealType mean, RealType p, RealType* result, const Policy& pol) |
| { |
| if((check_dist(function, mean, result, pol) == false) || |
| (check_prob(function, p, result, pol) == false)) |
| { |
| return false; |
| } |
| return true; |
| } // bool check_dist_and_prob |
| |
| } // namespace poisson_detail |
| |
| template <class RealType = double, class Policy = policies::policy<> > |
| class poisson_distribution |
| { |
| public: |
| typedef RealType value_type; |
| typedef Policy policy_type; |
| |
| poisson_distribution(RealType mean = 1) : m_l(mean) // mean (lambda). |
| { // Expected mean number of events that occur during the given interval. |
| RealType r; |
| poisson_detail::check_dist( |
| "boost::math::poisson_distribution<%1%>::poisson_distribution", |
| m_l, |
| &r, Policy()); |
| } // poisson_distribution constructor. |
| |
| RealType mean() const |
| { // Private data getter function. |
| return m_l; |
| } |
| private: |
| // Data member, initialized by constructor. |
| RealType m_l; // mean number of occurrences. |
| }; // template <class RealType, class Policy> class poisson_distribution |
| |
| typedef poisson_distribution<double> poisson; // Reserved name of type double. |
| |
| // Non-member functions to give properties of the distribution. |
| |
| template <class RealType, class Policy> |
| inline const std::pair<RealType, RealType> range(const poisson_distribution<RealType, Policy>& /* dist */) |
| { // Range of permissible values for random variable k. |
| using boost::math::tools::max_value; |
| return std::pair<RealType, RealType>(0, max_value<RealType>()); // Max integer? |
| } |
| |
| template <class RealType, class Policy> |
| inline const std::pair<RealType, RealType> support(const poisson_distribution<RealType, Policy>& /* dist */) |
| { // Range of supported values for random variable k. |
| // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. |
| using boost::math::tools::max_value; |
| return std::pair<RealType, RealType>(0, max_value<RealType>()); |
| } |
| |
| template <class RealType, class Policy> |
| inline RealType mean(const poisson_distribution<RealType, Policy>& dist) |
| { // Mean of poisson distribution = lambda. |
| return dist.mean(); |
| } // mean |
| |
| template <class RealType, class Policy> |
| inline RealType mode(const poisson_distribution<RealType, Policy>& dist) |
| { // mode. |
| BOOST_MATH_STD_USING // ADL of std functions. |
| return floor(dist.mean()); |
| } |
| |
| //template <class RealType, class Policy> |
| //inline RealType median(const poisson_distribution<RealType, Policy>& dist) |
| //{ // median = approximately lambda + 1/3 - 0.2/lambda |
| // RealType l = dist.mean(); |
| // return dist.mean() + static_cast<RealType>(0.3333333333333333333333333333333333333333333333) |
| // - static_cast<RealType>(0.2) / l; |
| //} // BUT this formula appears to be out-by-one compared to quantile(half) |
| // Query posted on Wikipedia. |
| // Now implemented via quantile(half) in derived accessors. |
| |
| template <class RealType, class Policy> |
| inline RealType variance(const poisson_distribution<RealType, Policy>& dist) |
| { // variance. |
| return dist.mean(); |
| } |
| |
| // RealType standard_deviation(const poisson_distribution<RealType, Policy>& dist) |
| // standard_deviation provided by derived accessors. |
| |
| template <class RealType, class Policy> |
| inline RealType skewness(const poisson_distribution<RealType, Policy>& dist) |
| { // skewness = sqrt(l). |
| BOOST_MATH_STD_USING // ADL of std functions. |
| return 1 / sqrt(dist.mean()); |
| } |
| |
| template <class RealType, class Policy> |
| inline RealType kurtosis_excess(const poisson_distribution<RealType, Policy>& dist) |
| { // skewness = sqrt(l). |
| return 1 / dist.mean(); // kurtosis_excess 1/mean from Wiki & MathWorld eq 31. |
| // http://mathworld.wolfram.com/Kurtosis.html explains that the kurtosis excess |
| // is more convenient because the kurtosis excess of a normal distribution is zero |
| // whereas the true kurtosis is 3. |
| } // RealType kurtosis_excess |
| |
| template <class RealType, class Policy> |
| inline RealType kurtosis(const poisson_distribution<RealType, Policy>& dist) |
| { // kurtosis is 4th moment about the mean = u4 / sd ^ 4 |
| // http://en.wikipedia.org/wiki/Curtosis |
| // kurtosis can range from -2 (flat top) to +infinity (sharp peak & heavy tails). |
| // http://www.itl.nist.gov/div898/handbook/eda/section3/eda35b.htm |
| return 3 + 1 / dist.mean(); // NIST. |
| // http://mathworld.wolfram.com/Kurtosis.html explains that the kurtosis excess |
| // is more convenient because the kurtosis excess of a normal distribution is zero |
| // whereas the true kurtosis is 3. |
| } // RealType kurtosis |
| |
| template <class RealType, class Policy> |
| RealType pdf(const poisson_distribution<RealType, Policy>& dist, const RealType& k) |
| { // Probability Density/Mass Function. |
| // Probability that there are EXACTLY k occurrences (or arrivals). |
| BOOST_FPU_EXCEPTION_GUARD |
| |
| BOOST_MATH_STD_USING // for ADL of std functions. |
| |
| RealType mean = dist.mean(); |
| // Error check: |
| RealType result; |
| if(false == poisson_detail::check_dist_and_k( |
| "boost::math::pdf(const poisson_distribution<%1%>&, %1%)", |
| mean, |
| k, |
| &result, Policy())) |
| { |
| return result; |
| } |
| |
| // Special case of mean zero, regardless of the number of events k. |
| if (mean == 0) |
| { // Probability for any k is zero. |
| return 0; |
| } |
| if (k == 0) |
| { // mean ^ k = 1, and k! = 1, so can simplify. |
| return exp(-mean); |
| } |
| return boost::math::gamma_p_derivative(k+1, mean, Policy()); |
| } // pdf |
| |
| template <class RealType, class Policy> |
| RealType cdf(const poisson_distribution<RealType, Policy>& dist, const RealType& k) |
| { // Cumulative Distribution Function Poisson. |
| // The random variate k is the number of occurrences(or arrivals) |
| // k argument may be integral, signed, or unsigned, or floating point. |
| // If necessary, it has already been promoted from an integral type. |
| // Returns the sum of the terms 0 through k of the Poisson Probability Density or Mass (pdf). |
| |
| // But note that the Poisson distribution |
| // (like others including the binomial, negative binomial & Bernoulli) |
| // is strictly defined as a discrete function: only integral values of k are envisaged. |
| // However because of the method of calculation using a continuous gamma function, |
| // it is convenient to treat it as if it is a continous function |
| // and permit non-integral values of k. |
| // To enforce the strict mathematical model, users should use floor or ceil functions |
| // outside this function to ensure that k is integral. |
| |
| // The terms are not summed directly (at least for larger k) |
| // instead the incomplete gamma integral is employed, |
| |
| BOOST_MATH_STD_USING // for ADL of std function exp. |
| |
| RealType mean = dist.mean(); |
| // Error checks: |
| RealType result; |
| if(false == poisson_detail::check_dist_and_k( |
| "boost::math::cdf(const poisson_distribution<%1%>&, %1%)", |
| mean, |
| k, |
| &result, Policy())) |
| { |
| return result; |
| } |
| // Special cases: |
| if (mean == 0) |
| { // Probability for any k is zero. |
| return 0; |
| } |
| if (k == 0) |
| { // return pdf(dist, static_cast<RealType>(0)); |
| // but mean (and k) have already been checked, |
| // so this avoids unnecessary repeated checks. |
| return exp(-mean); |
| } |
| // For small integral k could use a finite sum - |
| // it's cheaper than the gamma function. |
| // BUT this is now done efficiently by gamma_q function. |
| // Calculate poisson cdf using the gamma_q function. |
| return gamma_q(k+1, mean, Policy()); |
| } // binomial cdf |
| |
| template <class RealType, class Policy> |
| RealType cdf(const complemented2_type<poisson_distribution<RealType, Policy>, RealType>& c) |
| { // Complemented Cumulative Distribution Function Poisson |
| // The random variate k is the number of events, occurrences or arrivals. |
| // k argument may be integral, signed, or unsigned, or floating point. |
| // If necessary, it has already been promoted from an integral type. |
| // But note that the Poisson distribution |
| // (like others including the binomial, negative binomial & Bernoulli) |
| // is strictly defined as a discrete function: only integral values of k are envisaged. |
| // However because of the method of calculation using a continuous gamma function, |
| // it is convenient to treat it as is it is a continous function |
| // and permit non-integral values of k. |
| // To enforce the strict mathematical model, users should use floor or ceil functions |
| // outside this function to ensure that k is integral. |
| |
| // Returns the sum of the terms k+1 through inf of the Poisson Probability Density/Mass (pdf). |
| // The terms are not summed directly (at least for larger k) |
| // instead the incomplete gamma integral is employed, |
| |
| RealType const& k = c.param; |
| poisson_distribution<RealType, Policy> const& dist = c.dist; |
| |
| RealType mean = dist.mean(); |
| |
| // Error checks: |
| RealType result; |
| if(false == poisson_detail::check_dist_and_k( |
| "boost::math::cdf(const poisson_distribution<%1%>&, %1%)", |
| mean, |
| k, |
| &result, Policy())) |
| { |
| return result; |
| } |
| // Special case of mean, regardless of the number of events k. |
| if (mean == 0) |
| { // Probability for any k is unity, complement of zero. |
| return 1; |
| } |
| if (k == 0) |
| { // Avoid repeated checks on k and mean in gamma_p. |
| return -boost::math::expm1(-mean, Policy()); |
| } |
| // Unlike un-complemented cdf (sum from 0 to k), |
| // can't use finite sum from k+1 to infinity for small integral k, |
| // anyway it is now done efficiently by gamma_p. |
| return gamma_p(k + 1, mean, Policy()); // Calculate Poisson cdf using the gamma_p function. |
| // CCDF = gamma_p(k+1, lambda) |
| } // poisson ccdf |
| |
| template <class RealType, class Policy> |
| inline RealType quantile(const poisson_distribution<RealType, Policy>& dist, const RealType& p) |
| { // Quantile (or Percent Point) Poisson function. |
| // Return the number of expected events k for a given probability p. |
| RealType result; // of Argument checks: |
| if(false == poisson_detail::check_prob( |
| "boost::math::quantile(const poisson_distribution<%1%>&, %1%)", |
| p, |
| &result, Policy())) |
| { |
| return result; |
| } |
| // Special case: |
| if (dist.mean() == 0) |
| { // if mean = 0 then p = 0, so k can be anything? |
| if (false == poisson_detail::check_mean_NZ( |
| "boost::math::quantile(const poisson_distribution<%1%>&, %1%)", |
| dist.mean(), |
| &result, Policy())) |
| { |
| return result; |
| } |
| } |
| /* |
| BOOST_MATH_STD_USING // ADL of std functions. |
| // if(p == 0) NOT necessarily zero! |
| // Not necessarily any special value of k because is unlimited. |
| if (p <= exp(-dist.mean())) |
| { // if p <= cdf for 0 events (== pdf for 0 events), then quantile must be zero. |
| return 0; |
| } |
| return gamma_q_inva(dist.mean(), p, Policy()) - 1; |
| */ |
| typedef typename Policy::discrete_quantile_type discrete_type; |
| boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); |
| RealType guess, factor = 8; |
| RealType z = dist.mean(); |
| if(z < 1) |
| guess = z; |
| else |
| guess = boost::math::detail::inverse_poisson_cornish_fisher(z, p, RealType(1-p), Policy()); |
| if(z > 5) |
| { |
| if(z > 1000) |
| factor = 1.01f; |
| else if(z > 50) |
| factor = 1.1f; |
| else if(guess > 10) |
| factor = 1.25f; |
| else |
| factor = 2; |
| if(guess < 1.1) |
| factor = 8; |
| } |
| |
| return detail::inverse_discrete_quantile( |
| dist, |
| p, |
| 1-p, |
| guess, |
| factor, |
| RealType(1), |
| discrete_type(), |
| max_iter); |
| } // quantile |
| |
| template <class RealType, class Policy> |
| inline RealType quantile(const complemented2_type<poisson_distribution<RealType, Policy>, RealType>& c) |
| { // Quantile (or Percent Point) of Poisson function. |
| // Return the number of expected events k for a given |
| // complement of the probability q. |
| // |
| // Error checks: |
| RealType q = c.param; |
| const poisson_distribution<RealType, Policy>& dist = c.dist; |
| RealType result; // of argument checks. |
| if(false == poisson_detail::check_prob( |
| "boost::math::quantile(const poisson_distribution<%1%>&, %1%)", |
| q, |
| &result, Policy())) |
| { |
| return result; |
| } |
| // Special case: |
| if (dist.mean() == 0) |
| { // if mean = 0 then p = 0, so k can be anything? |
| if (false == poisson_detail::check_mean_NZ( |
| "boost::math::quantile(const poisson_distribution<%1%>&, %1%)", |
| dist.mean(), |
| &result, Policy())) |
| { |
| return result; |
| } |
| } |
| /* |
| if (-q <= boost::math::expm1(-dist.mean())) |
| { // if q <= cdf(complement for 0 events, then quantile must be zero. |
| return 0; |
| } |
| return gamma_p_inva(dist.mean(), q, Policy()) -1; |
| */ |
| typedef typename Policy::discrete_quantile_type discrete_type; |
| boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); |
| RealType guess, factor = 8; |
| RealType z = dist.mean(); |
| if(z < 1) |
| guess = z; |
| else |
| guess = boost::math::detail::inverse_poisson_cornish_fisher(z, RealType(1-q), q, Policy()); |
| if(z > 5) |
| { |
| if(z > 1000) |
| factor = 1.01f; |
| else if(z > 50) |
| factor = 1.1f; |
| else if(guess > 10) |
| factor = 1.25f; |
| else |
| factor = 2; |
| if(guess < 1.1) |
| factor = 8; |
| } |
| |
| return detail::inverse_discrete_quantile( |
| dist, |
| 1-q, |
| q, |
| guess, |
| factor, |
| RealType(1), |
| discrete_type(), |
| max_iter); |
| } // quantile complement. |
| |
| } // namespace math |
| } // namespace boost |
| |
| // This include must be at the end, *after* the accessors |
| // for this distribution have been defined, in order to |
| // keep compilers that support two-phase lookup happy. |
| #include <boost/math/distributions/detail/derived_accessors.hpp> |
| #include <boost/math/distributions/detail/inv_discrete_quantile.hpp> |
| |
| #endif // BOOST_MATH_SPECIAL_POISSON_HPP |
| |
| |
| |