572 lines · plain
1// boost\math\distributions\geometric.hpp2 3// Copyright John Maddock 2010.4// Copyright Paul A. Bristow 2010.5 6// Use, modification and distribution are subject to the7// Boost Software License, Version 1.0.8// (See accompanying file LICENSE_1_0.txt9// or copy at http://www.boost.org/LICENSE_1_0.txt)10 11// geometric distribution is a discrete probability distribution.12// It expresses the probability distribution of the number (k) of13// events, occurrences, failures or arrivals before the first success.14// supported on the set {0, 1, 2, 3...}15 16// Note that the set includes zero (unlike some definitions that start at one).17 18// The random variate k is the number of events, occurrences or arrivals.19// k argument may be integral, signed, or unsigned, or floating point.20// If necessary, it has already been promoted from an integral type.21 22// Note that the geometric distribution23// (like others including the binomial, geometric & Bernoulli)24// is strictly defined as a discrete function:25// only integral values of k are envisaged.26// However because the method of calculation uses a continuous gamma function,27// it is convenient to treat it as if a continuous function,28// and permit non-integral values of k.29// To enforce the strict mathematical model, users should use floor or ceil functions30// on k outside this function to ensure that k is integral.31 32// See http://en.wikipedia.org/wiki/geometric_distribution33// http://documents.wolfram.com/v5/Add-onsLinks/StandardPackages/Statistics/DiscreteDistributions.html34// http://mathworld.wolfram.com/GeometricDistribution.html35 36#ifndef BOOST_MATH_SPECIAL_GEOMETRIC_HPP37#define BOOST_MATH_SPECIAL_GEOMETRIC_HPP38 39#include <boost/math/tools/config.hpp>40#include <boost/math/tools/tuple.hpp>41#include <boost/math/tools/numeric_limits.hpp>42#include <boost/math/distributions/fwd.hpp>43#include <boost/math/special_functions/beta.hpp> // for ibeta(a, b, x) == Ix(a, b).44#include <boost/math/distributions/complement.hpp> // complement.45#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks domain_error & logic_error.46#include <boost/math/special_functions/fpclassify.hpp> // isnan.47#include <boost/math/tools/roots.hpp> // for root finding.48#include <boost/math/distributions/detail/inv_discrete_quantile.hpp>49#include <boost/math/special_functions/log1p.hpp>50 51#if defined (BOOST_MSVC)52# pragma warning(push)53// This believed not now necessary, so commented out.54//# pragma warning(disable: 4702) // unreachable code.55// in domain_error_imp in error_handling.56#endif57 58namespace boost59{60 namespace math61 {62 namespace geometric_detail63 {64 // Common error checking routines for geometric distribution function:65 template <class RealType, class Policy>66 BOOST_MATH_GPU_ENABLED inline bool check_success_fraction(const char* function, const RealType& p, RealType* result, const Policy& pol)67 {68 if( !(boost::math::isfinite)(p) || (p < 0) || (p > 1) )69 {70 *result = policies::raise_domain_error<RealType>(71 function,72 "Success fraction argument is %1%, but must be >= 0 and <= 1 !", p, pol);73 return false;74 }75 return true;76 }77 78 template <class RealType, class Policy>79 BOOST_MATH_GPU_ENABLED inline bool check_dist(const char* function, const RealType& p, RealType* result, const Policy& pol)80 {81 return check_success_fraction(function, p, result, pol);82 }83 84 template <class RealType, class Policy>85 BOOST_MATH_GPU_ENABLED inline bool check_dist_and_k(const char* function, const RealType& p, RealType k, RealType* result, const Policy& pol)86 {87 if(check_dist(function, p, result, pol) == false)88 {89 return false;90 }91 if( !(boost::math::isfinite)(k) || (k < 0) )92 { // Check k failures.93 *result = policies::raise_domain_error<RealType>(94 function,95 "Number of failures argument is %1%, but must be >= 0 !", k, pol);96 return false;97 }98 return true;99 } // Check_dist_and_k100 101 template <class RealType, class Policy>102 BOOST_MATH_GPU_ENABLED inline bool check_dist_and_prob(const char* function, RealType p, RealType prob, RealType* result, const Policy& pol)103 {104 if((check_dist(function, p, result, pol) && detail::check_probability(function, prob, result, pol)) == false)105 {106 return false;107 }108 return true;109 } // check_dist_and_prob110 } // namespace geometric_detail111 112 template <class RealType = double, class Policy = policies::policy<> >113 class geometric_distribution114 {115 public:116 typedef RealType value_type;117 typedef Policy policy_type;118 119 BOOST_MATH_GPU_ENABLED geometric_distribution(RealType p) : m_p(p)120 { // Constructor stores success_fraction p.121 RealType result;122 geometric_detail::check_dist(123 "geometric_distribution<%1%>::geometric_distribution",124 m_p, // Check success_fraction 0 <= p <= 1.125 &result, Policy());126 } // geometric_distribution constructor.127 128 // Private data getter class member functions.129 BOOST_MATH_GPU_ENABLED RealType success_fraction() const130 { // Probability of success as fraction in range 0 to 1.131 return m_p;132 }133 BOOST_MATH_GPU_ENABLED RealType successes() const134 { // Total number of successes r = 1 (for compatibility with negative binomial?).135 return 1;136 }137 138 // Parameter estimation.139 // (These are copies of negative_binomial distribution with successes = 1).140 BOOST_MATH_GPU_ENABLED static RealType find_lower_bound_on_p(141 RealType trials,142 RealType alpha) // alpha 0.05 equivalent to 95% for one-sided test.143 {144 constexpr auto function = "boost::math::geometric<%1%>::find_lower_bound_on_p";145 RealType result = 0; // of error checks.146 RealType successes = 1;147 RealType failures = trials - successes;148 if(false == detail::check_probability(function, alpha, &result, Policy())149 && geometric_detail::check_dist_and_k(150 function, RealType(0), failures, &result, Policy()))151 {152 return result;153 }154 // Use complement ibeta_inv function for lower bound.155 // This is adapted from the corresponding binomial formula156 // here: http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm157 // This is a Clopper-Pearson interval, and may be overly conservative,158 // see also "A Simple Improved Inferential Method for Some159 // Discrete Distributions" Yong CAI and K. KRISHNAMOORTHY160 // http://www.ucs.louisiana.edu/~kxk4695/Discrete_new.pdf161 //162 return ibeta_inv(successes, failures + 1, alpha, static_cast<RealType*>(nullptr), Policy());163 } // find_lower_bound_on_p164 165 BOOST_MATH_GPU_ENABLED static RealType find_upper_bound_on_p(166 RealType trials,167 RealType alpha) // alpha 0.05 equivalent to 95% for one-sided test.168 {169 constexpr auto function = "boost::math::geometric<%1%>::find_upper_bound_on_p";170 RealType result = 0; // of error checks.171 RealType successes = 1;172 RealType failures = trials - successes;173 if(false == geometric_detail::check_dist_and_k(174 function, RealType(0), failures, &result, Policy())175 && detail::check_probability(function, alpha, &result, Policy()))176 {177 return result;178 }179 if(failures == 0)180 {181 return 1;182 }// Use complement ibetac_inv function for upper bound.183 // Note adjusted failures value: *not* failures+1 as usual.184 // This is adapted from the corresponding binomial formula185 // here: http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm186 // This is a Clopper-Pearson interval, and may be overly conservative,187 // see also "A Simple Improved Inferential Method for Some188 // Discrete Distributions" Yong CAI and K. Krishnamoorthy189 // http://www.ucs.louisiana.edu/~kxk4695/Discrete_new.pdf190 //191 return ibetac_inv(successes, failures, alpha, static_cast<RealType*>(nullptr), Policy());192 } // find_upper_bound_on_p193 194 // Estimate number of trials :195 // "How many trials do I need to be P% sure of seeing k or fewer failures?"196 197 BOOST_MATH_GPU_ENABLED static RealType find_minimum_number_of_trials(198 RealType k, // number of failures (k >= 0).199 RealType p, // success fraction 0 <= p <= 1.200 RealType alpha) // risk level threshold 0 <= alpha <= 1.201 {202 constexpr auto function = "boost::math::geometric<%1%>::find_minimum_number_of_trials";203 // Error checks:204 RealType result = 0;205 if(false == geometric_detail::check_dist_and_k(206 function, p, k, &result, Policy())207 && detail::check_probability(function, alpha, &result, Policy()))208 {209 return result;210 }211 result = ibeta_inva(k + 1, p, alpha, Policy()); // returns n - k212 return result + k;213 } // RealType find_number_of_failures214 215 BOOST_MATH_GPU_ENABLED static RealType find_maximum_number_of_trials(216 RealType k, // number of failures (k >= 0).217 RealType p, // success fraction 0 <= p <= 1.218 RealType alpha) // risk level threshold 0 <= alpha <= 1.219 {220 constexpr auto function = "boost::math::geometric<%1%>::find_maximum_number_of_trials";221 // Error checks:222 RealType result = 0;223 if(false == geometric_detail::check_dist_and_k(224 function, p, k, &result, Policy())225 && detail::check_probability(function, alpha, &result, Policy()))226 {227 return result;228 }229 result = ibetac_inva(k + 1, p, alpha, Policy()); // returns n - k230 return result + k;231 } // RealType find_number_of_trials complemented232 233 private:234 //RealType m_r; // successes fixed at unity.235 RealType m_p; // success_fraction236 }; // template <class RealType, class Policy> class geometric_distribution237 238 typedef geometric_distribution<double> geometric; // Reserved name of type double.239 240 #ifdef __cpp_deduction_guides241 template <class RealType>242 geometric_distribution(RealType)->geometric_distribution<typename boost::math::tools::promote_args<RealType>::type>;243 #endif244 245 template <class RealType, class Policy>246 BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> range(const geometric_distribution<RealType, Policy>& /* dist */)247 { // Range of permissible values for random variable k.248 using boost::math::tools::max_value;249 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // max_integer?250 }251 252 template <class RealType, class Policy>253 BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> support(const geometric_distribution<RealType, Policy>& /* dist */)254 { // Range of supported values for random variable k.255 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.256 using boost::math::tools::max_value;257 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // max_integer?258 }259 260 template <class RealType, class Policy>261 BOOST_MATH_GPU_ENABLED inline RealType mean(const geometric_distribution<RealType, Policy>& dist)262 { // Mean of geometric distribution = (1-p)/p.263 return (1 - dist.success_fraction() ) / dist.success_fraction();264 } // mean265 266 // median implemented via quantile(half) in derived accessors.267 268 template <class RealType, class Policy>269 BOOST_MATH_GPU_ENABLED inline RealType mode(const geometric_distribution<RealType, Policy>&)270 { // Mode of geometric distribution = zero.271 BOOST_MATH_STD_USING // ADL of std functions.272 return 0;273 } // mode274 275 template <class RealType, class Policy>276 BOOST_MATH_GPU_ENABLED inline RealType variance(const geometric_distribution<RealType, Policy>& dist)277 { // Variance of Binomial distribution = (1-p) / p^2.278 return (1 - dist.success_fraction())279 / (dist.success_fraction() * dist.success_fraction());280 } // variance281 282 template <class RealType, class Policy>283 BOOST_MATH_GPU_ENABLED inline RealType skewness(const geometric_distribution<RealType, Policy>& dist)284 { // skewness of geometric distribution = 2-p / (sqrt(r(1-p))285 BOOST_MATH_STD_USING // ADL of std functions.286 RealType p = dist.success_fraction();287 return (2 - p) / sqrt(1 - p);288 } // skewness289 290 template <class RealType, class Policy>291 BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const geometric_distribution<RealType, Policy>& dist)292 { // kurtosis of geometric distribution293 // http://en.wikipedia.org/wiki/geometric is kurtosis_excess so add 3294 RealType p = dist.success_fraction();295 return 3 + (p*p - 6*p + 6) / (1 - p);296 } // kurtosis297 298 template <class RealType, class Policy>299 BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const geometric_distribution<RealType, Policy>& dist)300 { // kurtosis excess of geometric distribution301 // http://mathworld.wolfram.com/Kurtosis.html table of kurtosis_excess302 RealType p = dist.success_fraction();303 return (p*p - 6*p + 6) / (1 - p);304 } // kurtosis_excess305 306 // RealType standard_deviation(const geometric_distribution<RealType, Policy>& dist)307 // standard_deviation provided by derived accessors.308 // RealType hazard(const geometric_distribution<RealType, Policy>& dist)309 // hazard of geometric distribution provided by derived accessors.310 // RealType chf(const geometric_distribution<RealType, Policy>& dist)311 // chf of geometric distribution provided by derived accessors.312 313 template <class RealType, class Policy>314 BOOST_MATH_GPU_ENABLED inline RealType pdf(const geometric_distribution<RealType, Policy>& dist, const RealType& k)315 { // Probability Density/Mass Function.316 BOOST_FPU_EXCEPTION_GUARD317 BOOST_MATH_STD_USING // For ADL of math functions.318 constexpr auto function = "boost::math::pdf(const geometric_distribution<%1%>&, %1%)";319 320 RealType p = dist.success_fraction();321 RealType result = 0;322 if(false == geometric_detail::check_dist_and_k(323 function,324 p,325 k,326 &result, Policy()))327 {328 return result;329 }330 if (k == 0)331 {332 return p; // success_fraction333 }334 RealType q = 1 - p; // Inaccurate for small p?335 // So try to avoid inaccuracy for large or small p.336 // but has little effect > last significant bit.337 //cout << "p * pow(q, k) " << result << endl; // seems best whatever p338 //cout << "exp(p * k * log1p(-p)) " << p * exp(k * log1p(-p)) << endl;339 //if (p < 0.5)340 //{341 // result = p * pow(q, k);342 //}343 //else344 //{345 // result = p * exp(k * log1p(-p));346 //}347 result = p * pow(q, k);348 return result;349 } // geometric_pdf350 351 template <class RealType, class Policy>352 BOOST_MATH_GPU_ENABLED inline RealType cdf(const geometric_distribution<RealType, Policy>& dist, const RealType& k)353 { // Cumulative Distribution Function of geometric.354 constexpr auto function = "boost::math::cdf(const geometric_distribution<%1%>&, %1%)";355 356 // k argument may be integral, signed, or unsigned, or floating point.357 // If necessary, it has already been promoted from an integral type.358 RealType p = dist.success_fraction();359 // Error check:360 RealType result = 0;361 if(false == geometric_detail::check_dist_and_k(362 function,363 p,364 k,365 &result, Policy()))366 {367 return result;368 }369 if(k == 0)370 {371 return p; // success_fraction372 }373 //RealType q = 1 - p; // Bad for small p374 //RealType probability = 1 - std::pow(q, k+1);375 376 RealType z = boost::math::log1p(-p, Policy()) * (k + 1);377 RealType probability = -boost::math::expm1(z, Policy());378 379 return probability;380 } // cdf Cumulative Distribution Function geometric.381 382 template <class RealType, class Policy>383 BOOST_MATH_GPU_ENABLED inline RealType logcdf(const geometric_distribution<RealType, Policy>& dist, const RealType& k)384 { // Cumulative Distribution Function of geometric.385 BOOST_MATH_STD_USING386 constexpr auto function = "boost::math::logcdf(const geometric_distribution<%1%>&, %1%)";387 388 // k argument may be integral, signed, or unsigned, or floating point.389 // If necessary, it has already been promoted from an integral type.390 RealType p = dist.success_fraction();391 // Error check:392 RealType result = 0;393 if(false == geometric_detail::check_dist_and_k(394 function,395 p,396 k,397 &result, Policy()))398 {399 return -boost::math::numeric_limits<RealType>::infinity();400 }401 if(k == 0)402 {403 return log(p); // success_fraction404 }405 //RealType q = 1 - p; // Bad for small p406 //RealType probability = 1 - std::pow(q, k+1);407 408 RealType z = boost::math::log1p(-p, Policy()) * (k + 1);409 return log1p(-exp(z), Policy());410 } // logcdf Cumulative Distribution Function geometric.411 412 template <class RealType, class Policy>413 BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<geometric_distribution<RealType, Policy>, RealType>& c)414 { // Complemented Cumulative Distribution Function geometric.415 BOOST_MATH_STD_USING416 constexpr auto function = "boost::math::cdf(const geometric_distribution<%1%>&, %1%)";417 // k argument may be integral, signed, or unsigned, or floating point.418 // If necessary, it has already been promoted from an integral type.419 RealType const& k = c.param;420 geometric_distribution<RealType, Policy> const& dist = c.dist;421 RealType p = dist.success_fraction();422 // Error check:423 RealType result = 0;424 if(false == geometric_detail::check_dist_and_k(425 function,426 p,427 k,428 &result, Policy()))429 {430 return result;431 }432 RealType z = boost::math::log1p(-p, Policy()) * (k+1);433 RealType probability = exp(z);434 return probability;435 } // cdf Complemented Cumulative Distribution Function geometric.436 437 template <class RealType, class Policy>438 BOOST_MATH_GPU_ENABLED inline RealType logcdf(const complemented2_type<geometric_distribution<RealType, Policy>, RealType>& c)439 { // Complemented Cumulative Distribution Function geometric.440 BOOST_MATH_STD_USING441 constexpr auto function = "boost::math::logcdf(const geometric_distribution<%1%>&, %1%)";442 // k argument may be integral, signed, or unsigned, or floating point.443 // If necessary, it has already been promoted from an integral type.444 RealType const& k = c.param;445 geometric_distribution<RealType, Policy> const& dist = c.dist;446 RealType p = dist.success_fraction();447 // Error check:448 RealType result = 0;449 if(false == geometric_detail::check_dist_and_k(450 function,451 p,452 k,453 &result, Policy()))454 {455 return -boost::math::numeric_limits<RealType>::infinity();456 }457 458 return boost::math::log1p(-p, Policy()) * (k+1);459 } // logcdf Complemented Cumulative Distribution Function geometric.460 461 template <class RealType, class Policy>462 BOOST_MATH_GPU_ENABLED inline RealType quantile(const geometric_distribution<RealType, Policy>& dist, const RealType& x)463 { // Quantile, percentile/100 or Percent Point geometric function.464 // Return the number of expected failures k for a given probability p.465 466 // Inverse cumulative Distribution Function or Quantile (percentile / 100) of geometric Probability.467 // k argument may be integral, signed, or unsigned, or floating point.468 469 constexpr auto function = "boost::math::quantile(const geometric_distribution<%1%>&, %1%)";470 BOOST_MATH_STD_USING // ADL of std functions.471 472 RealType success_fraction = dist.success_fraction();473 // Check dist and x.474 RealType result = 0;475 if(false == geometric_detail::check_dist_and_prob476 (function, success_fraction, x, &result, Policy()))477 {478 return result;479 }480 481 // Special cases.482 if (x == 1)483 { // Would need +infinity failures for total confidence.484 result = policies::raise_overflow_error<RealType>(485 function,486 "Probability argument is 1, which implies infinite failures !", Policy());487 return result;488 // usually means return +std::numeric_limits<RealType>::infinity();489 // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR490 }491 if (x == 0)492 { // No failures are expected if P = 0.493 return 0; // Total trials will be just dist.successes.494 }495 // if (P <= pow(dist.success_fraction(), 1))496 if (x <= success_fraction)497 { // p <= pdf(dist, 0) == cdf(dist, 0)498 return 0;499 }500 if (x == 1)501 {502 return 0;503 }504 505 // log(1-x) /log(1-success_fraction) -1; but use log1p in case success_fraction is small506 result = boost::math::log1p(-x, Policy()) / boost::math::log1p(-success_fraction, Policy()) - 1;507 // Subtract a few epsilons here too?508 // to make sure it doesn't slip over, so ceil would be one too many.509 return result;510 } // RealType quantile(const geometric_distribution dist, p)511 512 template <class RealType, class Policy>513 BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<geometric_distribution<RealType, Policy>, RealType>& c)514 { // Quantile or Percent Point Binomial function.515 // Return the number of expected failures k for a given516 // complement of the probability Q = 1 - P.517 constexpr auto function = "boost::math::quantile(const geometric_distribution<%1%>&, %1%)";518 BOOST_MATH_STD_USING519 // Error checks:520 RealType x = c.param;521 const geometric_distribution<RealType, Policy>& dist = c.dist;522 RealType success_fraction = dist.success_fraction();523 RealType result = 0;524 if(false == geometric_detail::check_dist_and_prob(525 function,526 success_fraction,527 x,528 &result, Policy()))529 {530 return result;531 }532 533 // Special cases:534 if(x == 1)535 { // There may actually be no answer to this question,536 // since the probability of zero failures may be non-zero,537 return 0; // but zero is the best we can do:538 }539 if (-x <= boost::math::powm1(dist.success_fraction(), dist.successes(), Policy()))540 { // q <= cdf(complement(dist, 0)) == pdf(dist, 0)541 return 0; //542 }543 if(x == 0)544 { // Probability 1 - Q == 1 so infinite failures to achieve certainty.545 // Would need +infinity failures for total confidence.546 result = policies::raise_overflow_error<RealType>(547 function,548 "Probability argument complement is 0, which implies infinite failures !", Policy());549 return result;550 // usually means return +std::numeric_limits<RealType>::infinity();551 // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR552 }553 // log(x) /log(1-success_fraction) -1; but use log1p in case success_fraction is small554 result = log(x) / boost::math::log1p(-success_fraction, Policy()) - 1;555 return result;556 557 } // quantile complement558 559 } // namespace math560} // namespace boost561 562// This include must be at the end, *after* the accessors563// for this distribution have been defined, in order to564// keep compilers that support two-phase lookup happy.565#include <boost/math/distributions/detail/derived_accessors.hpp>566 567#if defined (BOOST_MSVC)568# pragma warning(pop)569#endif570 571#endif // BOOST_MATH_SPECIAL_GEOMETRIC_HPP572