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