730 lines · plain
1// boost\math\distributions\binomial.hpp2 3// Copyright John Maddock 2006.4// Copyright Paul A. Bristow 2007.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// http://en.wikipedia.org/wiki/binomial_distribution12 13// Binomial distribution is the discrete probability distribution of14// the number (k) of successes, in a sequence of15// n independent (yes or no, success or failure) Bernoulli trials.16 17// It expresses the probability of a number of events occurring in a fixed time18// if these events occur with a known average rate (probability of success),19// and are independent of the time since the last event.20 21// The number of cars that pass through a certain point on a road during a given period of time.22// The number of spelling mistakes a secretary makes while typing a single page.23// The number of phone calls at a call center per minute.24// The number of times a web server is accessed per minute.25// The number of light bulbs that burn out in a certain amount of time.26// The number of roadkill found per unit length of road27 28// http://en.wikipedia.org/wiki/binomial_distribution29 30// Given a sample of N measured values k[i],31// we wish to estimate the value of the parameter x (mean)32// of the binomial population from which the sample was drawn.33// To calculate the maximum likelihood value = 1/N sum i = 1 to N of k[i]34 35// Also may want a function for EXACTLY k.36 37// And probability that there are EXACTLY k occurrences is38// exp(-x) * pow(x, k) / factorial(k)39// where x is expected occurrences (mean) during the given interval.40// For example, if events occur, on average, every 4 min,41// and we are interested in number of events occurring in 10 min,42// then x = 10/4 = 2.543 44// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366i.htm45 46// The binomial distribution is used when there are47// exactly two mutually exclusive outcomes of a trial.48// These outcomes are appropriately labeled "success" and "failure".49// The binomial distribution is used to obtain50// the probability of observing x successes in N trials,51// with the probability of success on a single trial denoted by p.52// The binomial distribution assumes that p is fixed for all trials.53 54// P(x, p, n) = n!/(x! * (n-x)!) * p^x * (1-p)^(n-x)55 56// http://mathworld.wolfram.com/BinomialCoefficient.html57 58// The binomial coefficient (n; k) is the number of ways of picking59// k unordered outcomes from n possibilities,60// also known as a combination or combinatorial number.61// The symbols _nC_k and (n; k) are used to denote a binomial coefficient,62// and are sometimes read as "n choose k."63// (n; k) therefore gives the number of k-subsets possible out of a set of n distinct items.64 65// For example:66// The 2-subsets of {1,2,3,4} are the six pairs {1,2}, {1,3}, {1,4}, {2,3}, {2,4}, and {3,4}, so (4; 2)==6.67 68// http://functions.wolfram.com/GammaBetaErf/Binomial/ for evaluation.69 70// But note that the binomial distribution71// (like others including the poisson, negative binomial & Bernoulli)72// is strictly defined as a discrete function: only integral values of k are envisaged.73// However because of the method of calculation using a continuous gamma function,74// it is convenient to treat it as if a continuous function,75// and permit non-integral values of k.76// To enforce the strict mathematical model, users should use floor or ceil functions77// on k outside this function to ensure that k is integral.78 79#ifndef BOOST_MATH_SPECIAL_BINOMIAL_HPP80#define BOOST_MATH_SPECIAL_BINOMIAL_HPP81 82#include <boost/math/tools/config.hpp>83#include <boost/math/tools/tuple.hpp>84#include <boost/math/distributions/fwd.hpp>85#include <boost/math/special_functions/beta.hpp> // for incomplete beta.86#include <boost/math/distributions/complement.hpp> // complements87#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks88#include <boost/math/distributions/detail/inv_discrete_quantile.hpp> // error checks89#include <boost/math/special_functions/fpclassify.hpp> // isnan.90#include <boost/math/tools/roots.hpp> // for root finding.91 92#include <utility>93 94namespace boost95{96 namespace math97 {98 99 template <class RealType, class Policy>100 class binomial_distribution;101 102 namespace binomial_detail{103 // common error checking routines for binomial distribution functions:104 template <class RealType, class Policy>105 BOOST_MATH_CUDA_ENABLED inline bool check_N(const char* function, const RealType& N, RealType* result, const Policy& pol)106 {107 if((N < 0) || !(boost::math::isfinite)(N))108 {109 *result = policies::raise_domain_error<RealType>(110 function,111 "Number of Trials argument is %1%, but must be >= 0 !", N, pol);112 return false;113 }114 return true;115 }116 template <class RealType, class Policy>117 BOOST_MATH_CUDA_ENABLED inline bool check_success_fraction(const char* function, const RealType& p, RealType* result, const Policy& pol)118 {119 if((p < 0) || (p > 1) || !(boost::math::isfinite)(p))120 {121 *result = policies::raise_domain_error<RealType>(122 function,123 "Success fraction argument is %1%, but must be >= 0 and <= 1 !", p, pol);124 return false;125 }126 return true;127 }128 template <class RealType, class Policy>129 BOOST_MATH_CUDA_ENABLED inline bool check_dist(const char* function, const RealType& N, const RealType& p, RealType* result, const Policy& pol)130 {131 return check_success_fraction(132 function, p, result, pol)133 && check_N(134 function, N, result, pol);135 }136 template <class RealType, class Policy>137 BOOST_MATH_CUDA_ENABLED inline bool check_dist_and_k(const char* function, const RealType& N, const RealType& p, RealType k, RealType* result, const Policy& pol)138 {139 if(check_dist(function, N, p, result, pol) == false)140 return false;141 if((k < 0) || !(boost::math::isfinite)(k))142 {143 *result = policies::raise_domain_error<RealType>(144 function,145 "Number of Successes argument is %1%, but must be >= 0 !", k, pol);146 return false;147 }148 if(k > N)149 {150 *result = policies::raise_domain_error<RealType>(151 function,152 "Number of Successes argument is %1%, but must be <= Number of Trials !", k, pol);153 return false;154 }155 return true;156 }157 template <class RealType, class Policy>158 BOOST_MATH_CUDA_ENABLED inline bool check_dist_and_prob(const char* function, const RealType& N, RealType p, RealType prob, RealType* result, const Policy& pol)159 {160 if((check_dist(function, N, p, result, pol) && detail::check_probability(function, prob, result, pol)) == false)161 return false;162 return true;163 }164 165 template <class T, class Policy>166 BOOST_MATH_CUDA_ENABLED T inverse_binomial_cornish_fisher(T n, T sf, T p, T q, const Policy& pol)167 {168 BOOST_MATH_STD_USING169 // mean:170 T m = n * sf;171 // standard deviation:172 T sigma = sqrt(n * sf * (1 - sf));173 // skewness174 T sk = (1 - 2 * sf) / sigma;175 // kurtosis:176 // T k = (1 - 6 * sf * (1 - sf) ) / (n * sf * (1 - sf));177 // Get the inverse of a std normal distribution:178 T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();179 // Set the sign:180 if(p < 0.5)181 x = -x;182 T x2 = x * x;183 // w is correction term due to skewness184 T w = x + sk * (x2 - 1) / 6;185 /*186 // Add on correction due to kurtosis.187 // Disabled for now, seems to make things worse?188 //189 if(n >= 10)190 w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;191 */192 w = m + sigma * w;193 if(w < tools::min_value<T>())194 return sqrt(tools::min_value<T>());195 if(w > n)196 return n;197 return w;198 }199 200 template <class RealType, class Policy>201 BOOST_MATH_CUDA_ENABLED RealType quantile_imp(const binomial_distribution<RealType, Policy>& dist, const RealType& p, const RealType& q, bool comp)202 { // Quantile or Percent Point Binomial function.203 // Return the number of expected successes k,204 // for a given probability p.205 //206 // Error checks:207 BOOST_MATH_STD_USING // ADL of std names208 RealType result = 0;209 RealType trials = dist.trials();210 RealType success_fraction = dist.success_fraction();211 if(false == binomial_detail::check_dist_and_prob(212 "boost::math::quantile(binomial_distribution<%1%> const&, %1%)",213 trials,214 success_fraction,215 p,216 &result, Policy()))217 {218 return result;219 }220 221 // Special cases:222 //223 if(p == 0)224 { // There may actually be no answer to this question,225 // since the probability of zero successes may be non-zero,226 // but zero is the best we can do:227 return 0;228 }229 if(p == 1 || success_fraction == 1)230 { // Probability of n or fewer successes is always one,231 // so n is the most sensible answer here:232 return trials;233 }234 if (p <= pow(1 - success_fraction, trials))235 { // p <= pdf(dist, 0) == cdf(dist, 0)236 return 0; // So the only reasonable result is zero.237 } // And root finder would fail otherwise.238 239 // Solve for quantile numerically:240 //241 RealType guess = binomial_detail::inverse_binomial_cornish_fisher(trials, success_fraction, p, q, Policy());242 RealType factor = 8;243 if(trials > 100)244 factor = 1.01f; // guess is pretty accurate245 else if((trials > 10) && (trials - 1 > guess) && (guess > 3))246 factor = 1.15f; // less accurate but OK.247 else if(trials < 10)248 {249 // pretty inaccurate guess in this area:250 if(guess > trials / 64)251 {252 guess = trials / 4;253 factor = 2;254 }255 else256 guess = trials / 1024;257 }258 else259 factor = 2; // trials largish, but in far tails.260 261 typedef typename Policy::discrete_quantile_type discrete_quantile_type;262 std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();263 result = detail::inverse_discrete_quantile(264 dist,265 comp ? q : p,266 comp,267 guess,268 factor,269 RealType(1),270 discrete_quantile_type(),271 max_iter);272 return result;273 } // quantile274 275 }276 277 template <class RealType = double, class Policy = policies::policy<> >278 class binomial_distribution279 {280 public:281 typedef RealType value_type;282 typedef Policy policy_type;283 284 binomial_distribution(RealType n = 1, RealType p = 0.5) : m_n(n), m_p(p)285 { // Default n = 1 is the Bernoulli distribution286 // with equal probability of 'heads' or 'tails.287 RealType r;288 binomial_detail::check_dist(289 "boost::math::binomial_distribution<%1%>::binomial_distribution",290 m_n,291 m_p,292 &r, Policy());293 } // binomial_distribution constructor.294 295 BOOST_MATH_CUDA_ENABLED RealType success_fraction() const296 { // Probability.297 return m_p;298 }299 BOOST_MATH_CUDA_ENABLED RealType trials() const300 { // Total number of trials.301 return m_n;302 }303 304 enum interval_type{305 clopper_pearson_exact_interval,306 jeffreys_prior_interval307 };308 309 //310 // Estimation of the success fraction parameter.311 // The best estimate is actually simply successes/trials,312 // these functions are used313 // to obtain confidence intervals for the success fraction.314 //315 BOOST_MATH_CUDA_ENABLED static RealType find_lower_bound_on_p(316 RealType trials,317 RealType successes,318 RealType probability,319 interval_type t = clopper_pearson_exact_interval)320 {321 BOOST_MATH_STATIC const char* function = "boost::math::binomial_distribution<%1%>::find_lower_bound_on_p";322 // Error checks:323 RealType result = 0;324 if(false == binomial_detail::check_dist_and_k(325 function, trials, RealType(0), successes, &result, Policy())326 &&327 binomial_detail::check_dist_and_prob(328 function, trials, RealType(0), probability, &result, Policy()))329 { return result; }330 331 if(successes == 0)332 return 0;333 334 // NOTE!!! The Clopper Pearson formula uses "successes" not335 // "successes+1" as usual to get the lower bound,336 // see http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm337 return (t == clopper_pearson_exact_interval) ? ibeta_inv(successes, trials - successes + 1, probability, static_cast<RealType*>(nullptr), Policy())338 : ibeta_inv(successes + 0.5f, trials - successes + 0.5f, probability, static_cast<RealType*>(nullptr), Policy());339 }340 BOOST_MATH_CUDA_ENABLED static RealType find_upper_bound_on_p(341 RealType trials,342 RealType successes,343 RealType probability,344 interval_type t = clopper_pearson_exact_interval)345 {346 BOOST_MATH_STATIC const char* function = "boost::math::binomial_distribution<%1%>::find_upper_bound_on_p";347 // Error checks:348 RealType result = 0;349 if(false == binomial_detail::check_dist_and_k(350 function, trials, RealType(0), successes, &result, Policy())351 &&352 binomial_detail::check_dist_and_prob(353 function, trials, RealType(0), probability, &result, Policy()))354 { return result; }355 356 if(trials == successes)357 return 1;358 359 return (t == clopper_pearson_exact_interval) ? ibetac_inv(successes + 1, trials - successes, probability, static_cast<RealType*>(nullptr), Policy())360 : ibetac_inv(successes + 0.5f, trials - successes + 0.5f, probability, static_cast<RealType*>(nullptr), Policy());361 }362 // Estimate number of trials parameter:363 //364 // "How many trials do I need to be P% sure of seeing k events?"365 // or366 // "How many trials can I have to be P% sure of seeing fewer than k events?"367 //368 BOOST_MATH_CUDA_ENABLED static RealType find_minimum_number_of_trials(369 RealType k, // number of events370 RealType p, // success fraction371 RealType alpha) // risk level372 {373 BOOST_MATH_STATIC const char* function = "boost::math::binomial_distribution<%1%>::find_minimum_number_of_trials";374 // Error checks:375 RealType result = 0;376 if(false == binomial_detail::check_dist_and_k(377 function, k, p, k, &result, Policy())378 &&379 binomial_detail::check_dist_and_prob(380 function, k, p, alpha, &result, Policy()))381 { return result; }382 383 result = ibetac_invb(k + 1, p, alpha, Policy()); // returns n - k384 return result + k;385 }386 387 BOOST_MATH_CUDA_ENABLED static RealType find_maximum_number_of_trials(388 RealType k, // number of events389 RealType p, // success fraction390 RealType alpha) // risk level391 {392 BOOST_MATH_STATIC const char* function = "boost::math::binomial_distribution<%1%>::find_maximum_number_of_trials";393 // Error checks:394 RealType result = 0;395 if(false == binomial_detail::check_dist_and_k(396 function, k, p, k, &result, Policy())397 &&398 binomial_detail::check_dist_and_prob(399 function, k, p, alpha, &result, Policy()))400 { return result; }401 402 result = ibeta_invb(k + 1, p, alpha, Policy()); // returns n - k403 return result + k;404 }405 406 private:407 RealType m_n; // Not sure if this shouldn't be an int?408 RealType m_p; // success_fraction409 }; // template <class RealType, class Policy> class binomial_distribution410 411 typedef binomial_distribution<> binomial;412 // typedef binomial_distribution<double> binomial;413 // IS now included since no longer a name clash with function binomial.414 //typedef binomial_distribution<double> binomial; // Reserved name of type double.415 416 #ifdef __cpp_deduction_guides417 template <class RealType>418 binomial_distribution(RealType)->binomial_distribution<typename boost::math::tools::promote_args<RealType>::type>;419 template <class RealType>420 binomial_distribution(RealType,RealType)->binomial_distribution<typename boost::math::tools::promote_args<RealType>::type>;421 #endif422 423 template <class RealType, class Policy>424 BOOST_MATH_CUDA_ENABLED const boost::math::pair<RealType, RealType> range(const binomial_distribution<RealType, Policy>& dist)425 { // Range of permissible values for random variable k.426 using boost::math::tools::max_value;427 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), dist.trials());428 }429 430 template <class RealType, class Policy>431 BOOST_MATH_CUDA_ENABLED const boost::math::pair<RealType, RealType> support(const binomial_distribution<RealType, Policy>& dist)432 { // Range of supported values for random variable k.433 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.434 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), dist.trials());435 }436 437 template <class RealType, class Policy>438 BOOST_MATH_CUDA_ENABLED inline RealType mean(const binomial_distribution<RealType, Policy>& dist)439 { // Mean of Binomial distribution = np.440 return dist.trials() * dist.success_fraction();441 } // mean442 443 template <class RealType, class Policy>444 BOOST_MATH_CUDA_ENABLED inline RealType variance(const binomial_distribution<RealType, Policy>& dist)445 { // Variance of Binomial distribution = np(1-p).446 return dist.trials() * dist.success_fraction() * (1 - dist.success_fraction());447 } // variance448 449 template <class RealType, class Policy>450 BOOST_MATH_CUDA_ENABLED RealType pdf(const binomial_distribution<RealType, Policy>& dist, const RealType& k)451 { // Probability Density/Mass Function.452 BOOST_FPU_EXCEPTION_GUARD453 454 BOOST_MATH_STD_USING // for ADL of std functions455 456 RealType n = dist.trials();457 458 // Error check:459 RealType result = 0; // initialization silences some compiler warnings460 if(false == binomial_detail::check_dist_and_k(461 "boost::math::pdf(binomial_distribution<%1%> const&, %1%)",462 n,463 dist.success_fraction(),464 k,465 &result, Policy()))466 {467 return result;468 }469 470 // Special cases of success_fraction, regardless of k successes and regardless of n trials.471 if (dist.success_fraction() == 0)472 { // probability of zero successes is 1:473 return static_cast<RealType>(k == 0 ? 1 : 0);474 }475 if (dist.success_fraction() == 1)476 { // probability of n successes is 1:477 return static_cast<RealType>(k == n ? 1 : 0);478 }479 // k argument may be integral, signed, or unsigned, or floating point.480 // If necessary, it has already been promoted from an integral type.481 if (n == 0)482 {483 return 1; // Probability = 1 = certainty.484 }485 if (k == n)486 { // binomial coeffic (n n) = 1,487 // n ^ 0 = 1488 return pow(dist.success_fraction(), k); // * pow((1 - dist.success_fraction()), (n - k)) = 1489 }490 491 // Probability of getting exactly k successes492 // if C(n, k) is the binomial coefficient then:493 //494 // f(k; n,p) = C(n, k) * p^k * (1-p)^(n-k)495 // = (n!/(k!(n-k)!)) * p^k * (1-p)^(n-k)496 // = (tgamma(n+1) / (tgamma(k+1)*tgamma(n-k+1))) * p^k * (1-p)^(n-k)497 // = p^k (1-p)^(n-k) / (beta(k+1, n-k+1) * (n+1))498 // = ibeta_derivative(k+1, n-k+1, p) / (n+1)499 //500 using boost::math::ibeta_derivative; // a, b, x501 return ibeta_derivative(k+1, n-k+1, dist.success_fraction(), Policy()) / (n+1);502 503 } // pdf504 505 template <class RealType, class Policy>506 BOOST_MATH_CUDA_ENABLED inline RealType cdf(const binomial_distribution<RealType, Policy>& dist, const RealType& k)507 { // Cumulative Distribution Function Binomial.508 // The random variate k is the number of successes in n trials.509 // k argument may be integral, signed, or unsigned, or floating point.510 // If necessary, it has already been promoted from an integral type.511 512 // Returns the sum of the terms 0 through k of the Binomial Probability Density/Mass:513 //514 // i=k515 // -- ( n ) i n-i516 // > | | p (1-p)517 // -- ( i )518 // i=0519 520 // The terms are not summed directly instead521 // the incomplete beta integral is employed,522 // according to the formula:523 // P = I[1-p]( n-k, k+1).524 // = 1 - I[p](k + 1, n - k)525 526 BOOST_MATH_STD_USING // for ADL of std functions527 528 RealType n = dist.trials();529 RealType p = dist.success_fraction();530 531 // Error check:532 RealType result = 0;533 if(false == binomial_detail::check_dist_and_k(534 "boost::math::cdf(binomial_distribution<%1%> const&, %1%)",535 n,536 p,537 k,538 &result, Policy()))539 {540 return result;541 }542 if (k == n)543 {544 return 1;545 }546 547 // Special cases, regardless of k.548 if (p == 0)549 { // This need explanation:550 // the pdf is zero for all cases except when k == 0.551 // For zero p the probability of zero successes is one.552 // Therefore the cdf is always 1:553 // the probability of k or *fewer* successes is always 1554 // if there are never any successes!555 return 1;556 }557 if (p == 1)558 { // This is correct but needs explanation:559 // when k = 1560 // all the cdf and pdf values are zero *except* when k == n,561 // and that case has been handled above already.562 return 0;563 }564 //565 // P = I[1-p](n - k, k + 1)566 // = 1 - I[p](k + 1, n - k)567 // Use of ibetac here prevents cancellation errors in calculating568 // 1-p if p is very small, perhaps smaller than machine epsilon.569 //570 // Note that we do not use a finite sum here, since the incomplete571 // beta uses a finite sum internally for integer arguments, so572 // we'll just let it take care of the necessary logic.573 //574 return ibetac(k + 1, n - k, p, Policy());575 } // binomial cdf576 577 template <class RealType, class Policy>578 BOOST_MATH_CUDA_ENABLED inline RealType cdf(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c)579 { // Complemented Cumulative Distribution Function Binomial.580 // The random variate k is the number of successes in n trials.581 // k argument may be integral, signed, or unsigned, or floating point.582 // If necessary, it has already been promoted from an integral type.583 584 // Returns the sum of the terms k+1 through n of the Binomial Probability Density/Mass:585 //586 // i=n587 // -- ( n ) i n-i588 // > | | p (1-p)589 // -- ( i )590 // i=k+1591 592 // The terms are not summed directly instead593 // the incomplete beta integral is employed,594 // according to the formula:595 // Q = 1 -I[1-p]( n-k, k+1).596 // = I[p](k + 1, n - k)597 598 BOOST_MATH_STD_USING // for ADL of std functions599 600 RealType const& k = c.param;601 binomial_distribution<RealType, Policy> const& dist = c.dist;602 RealType n = dist.trials();603 RealType p = dist.success_fraction();604 605 // Error checks:606 RealType result = 0;607 if(false == binomial_detail::check_dist_and_k(608 "boost::math::cdf(binomial_distribution<%1%> const&, %1%)",609 n,610 p,611 k,612 &result, Policy()))613 {614 return result;615 }616 617 if (k == n)618 { // Probability of greater than n successes is necessarily zero:619 return 0;620 }621 622 // Special cases, regardless of k.623 if (p == 0)624 {625 // This need explanation: the pdf is zero for all626 // cases except when k == 0. For zero p the probability627 // of zero successes is one. Therefore the cdf is always628 // 1: the probability of *more than* k successes is always 0629 // if there are never any successes!630 return 0;631 }632 if (p == 1)633 {634 // This needs explanation, when p = 1635 // we always have n successes, so the probability636 // of more than k successes is 1 as long as k < n.637 // The k == n case has already been handled above.638 return 1;639 }640 //641 // Calculate cdf binomial using the incomplete beta function.642 // Q = 1 -I[1-p](n - k, k + 1)643 // = I[p](k + 1, n - k)644 // Use of ibeta here prevents cancellation errors in calculating645 // 1-p if p is very small, perhaps smaller than machine epsilon.646 //647 // Note that we do not use a finite sum here, since the incomplete648 // beta uses a finite sum internally for integer arguments, so649 // we'll just let it take care of the necessary logic.650 //651 return ibeta(k + 1, n - k, p, Policy());652 } // binomial cdf653 654 template <class RealType, class Policy>655 BOOST_MATH_CUDA_ENABLED inline RealType quantile(const binomial_distribution<RealType, Policy>& dist, const RealType& p)656 {657 return binomial_detail::quantile_imp(dist, p, RealType(1-p), false);658 } // quantile659 660 template <class RealType, class Policy>661 BOOST_MATH_CUDA_ENABLED RealType quantile(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c)662 {663 return binomial_detail::quantile_imp(c.dist, RealType(1-c.param), c.param, true);664 } // quantile665 666 template <class RealType, class Policy>667 BOOST_MATH_CUDA_ENABLED inline RealType mode(const binomial_distribution<RealType, Policy>& dist)668 {669 BOOST_MATH_STD_USING // ADL of std functions.670 RealType p = dist.success_fraction();671 RealType n = dist.trials();672 return floor(p * (n + 1));673 }674 675 template <class RealType, class Policy>676 BOOST_MATH_CUDA_ENABLED inline RealType median(const binomial_distribution<RealType, Policy>& dist)677 { // Bounds for the median of the negative binomial distribution678 // VAN DE VEN R. ; WEBER N. C. ;679 // Univ. Sydney, school mathematics statistics, Sydney N.S.W. 2006, AUSTRALIE680 // Metrika (Metrika) ISSN 0026-1335 CODEN MTRKA8681 // 1993, vol. 40, no3-4, pp. 185-189 (4 ref.)682 683 // Bounds for median and 50 percentage point of binomial and negative binomial distribution684 // Metrika, ISSN 0026-1335 (Print) 1435-926X (Online)685 // Volume 41, Number 1 / December, 1994, DOI 10.1007/BF01895303686 BOOST_MATH_STD_USING // ADL of std functions.687 RealType p = dist.success_fraction();688 RealType n = dist.trials();689 // Wikipedia says one of floor(np) -1, floor (np), floor(np) +1690 return floor(p * n); // Chose the middle value.691 }692 693 template <class RealType, class Policy>694 BOOST_MATH_CUDA_ENABLED inline RealType skewness(const binomial_distribution<RealType, Policy>& dist)695 {696 BOOST_MATH_STD_USING // ADL of std functions.697 RealType p = dist.success_fraction();698 RealType n = dist.trials();699 return (1 - 2 * p) / sqrt(n * p * (1 - p));700 }701 702 template <class RealType, class Policy>703 BOOST_MATH_CUDA_ENABLED inline RealType kurtosis(const binomial_distribution<RealType, Policy>& dist)704 {705 RealType p = dist.success_fraction();706 RealType n = dist.trials();707 return 3 - 6 / n + 1 / (n * p * (1 - p));708 }709 710 template <class RealType, class Policy>711 BOOST_MATH_CUDA_ENABLED inline RealType kurtosis_excess(const binomial_distribution<RealType, Policy>& dist)712 {713 RealType p = dist.success_fraction();714 RealType q = 1 - p;715 RealType n = dist.trials();716 return (1 - 6 * p * q) / (n * p * q);717 }718 719 } // namespace math720 } // namespace boost721 722// This include must be at the end, *after* the accessors723// for this distribution have been defined, in order to724// keep compilers that support two-phase lookup happy.725#include <boost/math/distributions/detail/derived_accessors.hpp>726 727#endif // BOOST_MATH_SPECIAL_BINOMIAL_HPP728 729 730