comparison any/include/boost/math/distributions/poisson.hpp @ 160:cff480c41f97

Add some cross-platform Boost headers

+// 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 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 l_mean = 1) : m_l(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>(static_cast<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>(static_cast<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 = 0;
+      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 = 0;
+      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 = 0;
+      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.
+      static const char* function = "boost::math::quantile(const poisson_distribution<%1%>&, %1%)";
+      RealType result = 0; // of Argument checks:
+      if(false == poisson_detail::check_prob(
+        function,
+        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(
+         function,
+         dist.mean(),
+         &result, Policy()))
+        {
+          return result;
+        }
+      }
+      if(p == 0)
+      {
+         return 0; // Exact result regardless of discrete-quantile Policy
+      }
+      if(p == 1)
+      {
+         return policies::raise_overflow_error<RealType>(function, 0, Policy());
+      }
+      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,
+         false,
+         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:
+      static const char* function = "boost::math::quantile(complement(const poisson_distribution<%1%>&, %1%))";
+      RealType q = c.param;
+      const poisson_distribution<RealType, Policy>& dist = c.dist;
+      RealType result = 0;  // of argument checks.
+      if(false == poisson_detail::check_prob(
+        function,
+        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(
+         function,
+         dist.mean(),
+         &result, Policy()))
+        {
+          return result;
+        }
+      }
+      if(q == 0)
+      {
+         return policies::raise_overflow_error<RealType>(function, 0, Policy());
+      }
+      if(q == 1)
+      {
+         return 0;  // Exact result regardless of discrete-quantile Policy
+      }
+      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,
+         q,
+         true,
+         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
+
+
+