Chris@16: // boost\math\distributions\binomial.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: // http://en.wikipedia.org/wiki/binomial_distribution Chris@16: Chris@16: // Binomial distribution is the discrete probability distribution of Chris@16: // the number (k) of successes, in a sequence of Chris@16: // n independent (yes or no, success or failure) Bernoulli trials. Chris@16: Chris@16: // It expresses the probability of a number of events occurring in a fixed time Chris@16: // if these events occur with a known average rate (probability of success), Chris@16: // and are independent of the time since the last event. Chris@16: Chris@16: // The number of cars that pass through a certain point on a road during a given period of time. Chris@16: // The number of spelling mistakes a secretary makes while typing a single page. Chris@16: // The number of phone calls at a call center per minute. Chris@16: // The number of times a web server is accessed per minute. Chris@16: // The number of light bulbs that burn out in a certain amount of time. Chris@16: // The number of roadkill found per unit length of road Chris@16: Chris@16: // http://en.wikipedia.org/wiki/binomial_distribution Chris@16: Chris@16: // Given a sample of N measured values k[i], Chris@16: // we wish to estimate the value of the parameter x (mean) Chris@16: // of the binomial population from which the sample was drawn. Chris@16: // To calculate the maximum likelihood value = 1/N sum i = 1 to N of k[i] Chris@16: Chris@16: // Also may want a function for EXACTLY k. Chris@16: Chris@16: // And probability that there are EXACTLY k occurrences is Chris@16: // exp(-x) * pow(x, k) / factorial(k) Chris@16: // where x is expected occurrences (mean) during the given interval. Chris@16: // For example, if events occur, on average, every 4 min, Chris@16: // and we are interested in number of events occurring in 10 min, Chris@16: // then x = 10/4 = 2.5 Chris@16: Chris@16: // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366i.htm Chris@16: Chris@16: // The binomial distribution is used when there are Chris@16: // exactly two mutually exclusive outcomes of a trial. Chris@16: // These outcomes are appropriately labeled "success" and "failure". Chris@16: // The binomial distribution is used to obtain Chris@16: // the probability of observing x successes in N trials, Chris@16: // with the probability of success on a single trial denoted by p. Chris@16: // The binomial distribution assumes that p is fixed for all trials. Chris@16: Chris@16: // P(x, p, n) = n!/(x! * (n-x)!) * p^x * (1-p)^(n-x) Chris@16: Chris@16: // http://mathworld.wolfram.com/BinomialCoefficient.html Chris@16: Chris@16: // The binomial coefficient (n; k) is the number of ways of picking Chris@16: // k unordered outcomes from n possibilities, Chris@16: // also known as a combination or combinatorial number. Chris@16: // The symbols _nC_k and (n; k) are used to denote a binomial coefficient, Chris@16: // and are sometimes read as "n choose k." Chris@16: // (n; k) therefore gives the number of k-subsets possible out of a set of n distinct items. Chris@16: Chris@16: // For example: Chris@16: // 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. Chris@16: Chris@16: // http://functions.wolfram.com/GammaBetaErf/Binomial/ for evaluation. Chris@16: Chris@16: // But note that the binomial distribution Chris@16: // (like others including the poisson, 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 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: #ifndef BOOST_MATH_SPECIAL_BINOMIAL_HPP Chris@16: #define BOOST_MATH_SPECIAL_BINOMIAL_HPP Chris@16: Chris@16: #include Chris@16: #include // for incomplete beta. Chris@16: #include // complements Chris@16: #include // error checks Chris@16: #include // error checks Chris@16: #include // isnan. Chris@16: #include // for root finding. Chris@16: Chris@16: #include Chris@16: Chris@16: namespace boost Chris@16: { Chris@16: namespace math Chris@16: { Chris@16: Chris@16: template Chris@16: class binomial_distribution; Chris@16: Chris@16: namespace binomial_detail{ Chris@16: // common error checking routines for binomial distribution functions: Chris@16: template Chris@16: inline bool check_N(const char* function, const RealType& N, RealType* result, const Policy& pol) Chris@16: { Chris@16: if((N < 0) || !(boost::math::isfinite)(N)) Chris@16: { Chris@16: *result = policies::raise_domain_error( Chris@16: function, Chris@16: "Number of Trials argument is %1%, but must be >= 0 !", N, pol); Chris@16: return false; Chris@16: } Chris@16: return true; Chris@16: } Chris@16: template Chris@16: inline bool check_success_fraction(const char* function, const RealType& p, RealType* result, const Policy& pol) Chris@16: { Chris@16: if((p < 0) || (p > 1) || !(boost::math::isfinite)(p)) Chris@16: { Chris@16: *result = policies::raise_domain_error( Chris@16: function, Chris@16: "Success fraction argument is %1%, but must be >= 0 and <= 1 !", p, pol); Chris@16: return false; Chris@16: } Chris@16: return true; Chris@16: } Chris@16: template Chris@16: inline bool check_dist(const char* function, const RealType& N, const RealType& p, RealType* result, const Policy& pol) Chris@16: { Chris@16: return check_success_fraction( Chris@16: function, p, result, pol) Chris@16: && check_N( Chris@16: function, N, result, pol); Chris@16: } Chris@16: template Chris@16: inline bool check_dist_and_k(const char* function, const RealType& N, const RealType& p, RealType k, RealType* result, const Policy& pol) Chris@16: { Chris@16: if(check_dist(function, N, p, result, pol) == false) Chris@16: return false; Chris@16: if((k < 0) || !(boost::math::isfinite)(k)) Chris@16: { Chris@16: *result = policies::raise_domain_error( Chris@16: function, Chris@16: "Number of Successes argument is %1%, but must be >= 0 !", k, pol); Chris@16: return false; Chris@16: } Chris@16: if(k > N) Chris@16: { Chris@16: *result = policies::raise_domain_error( Chris@16: function, Chris@16: "Number of Successes argument is %1%, but must be <= Number of Trials !", k, pol); Chris@16: return false; Chris@16: } Chris@16: return true; Chris@16: } Chris@16: template Chris@16: inline bool check_dist_and_prob(const char* function, const RealType& N, RealType p, RealType prob, RealType* result, const Policy& pol) Chris@16: { Chris@16: if(check_dist(function, N, p, result, pol) && detail::check_probability(function, prob, result, pol) == false) Chris@16: return false; Chris@16: return true; Chris@16: } Chris@16: Chris@16: template Chris@16: T inverse_binomial_cornish_fisher(T n, T sf, T p, T q, const Policy& pol) Chris@16: { Chris@16: BOOST_MATH_STD_USING Chris@16: // mean: Chris@16: T m = n * sf; Chris@16: // standard deviation: Chris@16: T sigma = sqrt(n * sf * (1 - sf)); Chris@16: // skewness Chris@16: T sk = (1 - 2 * sf) / sigma; Chris@16: // kurtosis: Chris@16: // T k = (1 - 6 * sf * (1 - sf) ) / (n * sf * (1 - sf)); Chris@16: // Get the inverse of a std normal distribution: Chris@16: T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two(); Chris@16: // Set the sign: Chris@16: if(p < 0.5) Chris@16: x = -x; Chris@16: T x2 = x * x; Chris@16: // w is correction term due to skewness Chris@16: T w = x + sk * (x2 - 1) / 6; Chris@16: /* Chris@16: // Add on correction due to kurtosis. Chris@16: // Disabled for now, seems to make things worse? Chris@16: // Chris@16: if(n >= 10) Chris@16: w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36; Chris@16: */ Chris@16: w = m + sigma * w; Chris@16: if(w < tools::min_value()) Chris@16: return sqrt(tools::min_value()); Chris@16: if(w > n) Chris@16: return n; Chris@16: return w; Chris@16: } Chris@16: Chris@16: template Chris@16: RealType quantile_imp(const binomial_distribution& dist, const RealType& p, const RealType& q, bool comp) Chris@16: { // Quantile or Percent Point Binomial function. Chris@16: // Return the number of expected successes k, Chris@16: // for a given probability p. Chris@16: // Chris@16: // Error checks: Chris@16: BOOST_MATH_STD_USING // ADL of std names Chris@16: RealType result = 0; Chris@16: RealType trials = dist.trials(); Chris@16: RealType success_fraction = dist.success_fraction(); Chris@16: if(false == binomial_detail::check_dist_and_prob( Chris@16: "boost::math::quantile(binomial_distribution<%1%> const&, %1%)", Chris@16: trials, Chris@16: success_fraction, Chris@16: p, Chris@16: &result, Policy())) Chris@16: { Chris@16: return result; Chris@16: } Chris@16: Chris@16: // Special cases: Chris@16: // Chris@16: if(p == 0) Chris@16: { // There may actually be no answer to this question, Chris@16: // since the probability of zero successes may be non-zero, Chris@16: // but zero is the best we can do: Chris@16: return 0; Chris@16: } Chris@16: if(p == 1) Chris@16: { // Probability of n or fewer successes is always one, Chris@16: // so n is the most sensible answer here: Chris@16: return trials; Chris@16: } Chris@16: if (p <= pow(1 - success_fraction, trials)) Chris@16: { // p <= pdf(dist, 0) == cdf(dist, 0) Chris@16: return 0; // So the only reasonable result is zero. Chris@16: } // And root finder would fail otherwise. Chris@16: if(success_fraction == 1) Chris@16: { // our formulae break down in this case: Chris@16: return p > 0.5f ? trials : 0; Chris@16: } Chris@16: Chris@16: // Solve for quantile numerically: Chris@16: // Chris@16: RealType guess = binomial_detail::inverse_binomial_cornish_fisher(trials, success_fraction, p, q, Policy()); Chris@16: RealType factor = 8; Chris@16: if(trials > 100) Chris@16: factor = 1.01f; // guess is pretty accurate Chris@16: else if((trials > 10) && (trials - 1 > guess) && (guess > 3)) Chris@16: factor = 1.15f; // less accurate but OK. Chris@16: else if(trials < 10) Chris@16: { Chris@16: // pretty inaccurate guess in this area: Chris@16: if(guess > trials / 64) Chris@16: { Chris@16: guess = trials / 4; Chris@16: factor = 2; Chris@16: } Chris@16: else Chris@16: guess = trials / 1024; Chris@16: } Chris@16: else Chris@16: factor = 2; // trials largish, but in far tails. Chris@16: Chris@16: typedef typename Policy::discrete_quantile_type discrete_quantile_type; Chris@16: boost::uintmax_t max_iter = policies::get_max_root_iterations(); Chris@16: return detail::inverse_discrete_quantile( Chris@16: dist, Chris@16: comp ? q : p, Chris@16: comp, Chris@16: guess, Chris@16: factor, Chris@16: RealType(1), Chris@16: discrete_quantile_type(), Chris@16: max_iter); Chris@16: } // quantile Chris@16: Chris@16: } Chris@16: Chris@16: template > Chris@16: class binomial_distribution Chris@16: { Chris@16: public: Chris@16: typedef RealType value_type; Chris@16: typedef Policy policy_type; Chris@16: Chris@16: binomial_distribution(RealType n = 1, RealType p = 0.5) : m_n(n), m_p(p) Chris@16: { // Default n = 1 is the Bernoulli distribution Chris@16: // with equal probability of 'heads' or 'tails. Chris@16: RealType r; Chris@16: binomial_detail::check_dist( Chris@16: "boost::math::binomial_distribution<%1%>::binomial_distribution", Chris@16: m_n, Chris@16: m_p, Chris@16: &r, Policy()); Chris@16: } // binomial_distribution constructor. Chris@16: Chris@16: RealType success_fraction() const Chris@16: { // Probability. Chris@16: return m_p; Chris@16: } Chris@16: RealType trials() const Chris@16: { // Total number of trials. Chris@16: return m_n; Chris@16: } Chris@16: Chris@16: enum interval_type{ Chris@16: clopper_pearson_exact_interval, Chris@16: jeffreys_prior_interval Chris@16: }; Chris@16: Chris@16: // Chris@16: // Estimation of the success fraction parameter. Chris@16: // The best estimate is actually simply successes/trials, Chris@16: // these functions are used Chris@16: // to obtain confidence intervals for the success fraction. Chris@16: // Chris@16: static RealType find_lower_bound_on_p( Chris@16: RealType trials, Chris@16: RealType successes, Chris@16: RealType probability, Chris@16: interval_type t = clopper_pearson_exact_interval) Chris@16: { Chris@16: static const char* function = "boost::math::binomial_distribution<%1%>::find_lower_bound_on_p"; Chris@16: // Error checks: Chris@16: RealType result = 0; Chris@16: if(false == binomial_detail::check_dist_and_k( Chris@16: function, trials, RealType(0), successes, &result, Policy()) Chris@16: && Chris@16: binomial_detail::check_dist_and_prob( Chris@16: function, trials, RealType(0), probability, &result, Policy())) Chris@16: { return result; } Chris@16: Chris@16: if(successes == 0) Chris@16: return 0; Chris@16: Chris@16: // NOTE!!! The Clopper Pearson formula uses "successes" not Chris@16: // "successes+1" as usual to get the lower bound, Chris@16: // see http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm Chris@16: return (t == clopper_pearson_exact_interval) ? ibeta_inv(successes, trials - successes + 1, probability, static_cast(0), Policy()) Chris@16: : ibeta_inv(successes + 0.5f, trials - successes + 0.5f, probability, static_cast(0), Policy()); Chris@16: } Chris@16: static RealType find_upper_bound_on_p( Chris@16: RealType trials, Chris@16: RealType successes, Chris@16: RealType probability, Chris@16: interval_type t = clopper_pearson_exact_interval) Chris@16: { Chris@16: static const char* function = "boost::math::binomial_distribution<%1%>::find_upper_bound_on_p"; Chris@16: // Error checks: Chris@16: RealType result = 0; Chris@16: if(false == binomial_detail::check_dist_and_k( Chris@16: function, trials, RealType(0), successes, &result, Policy()) Chris@16: && Chris@16: binomial_detail::check_dist_and_prob( Chris@16: function, trials, RealType(0), probability, &result, Policy())) Chris@16: { return result; } Chris@16: Chris@16: if(trials == successes) Chris@16: return 1; Chris@16: Chris@16: return (t == clopper_pearson_exact_interval) ? ibetac_inv(successes + 1, trials - successes, probability, static_cast(0), Policy()) Chris@16: : ibetac_inv(successes + 0.5f, trials - successes + 0.5f, probability, static_cast(0), Policy()); Chris@16: } Chris@16: // Estimate number of trials parameter: Chris@16: // Chris@16: // "How many trials do I need to be P% sure of seeing k events?" Chris@16: // or Chris@16: // "How many trials can I have to be P% sure of seeing fewer than k events?" Chris@16: // Chris@16: static RealType find_minimum_number_of_trials( Chris@16: RealType k, // number of events Chris@16: RealType p, // success fraction Chris@16: RealType alpha) // risk level Chris@16: { Chris@16: static const char* function = "boost::math::binomial_distribution<%1%>::find_minimum_number_of_trials"; Chris@16: // Error checks: Chris@16: RealType result = 0; Chris@16: if(false == binomial_detail::check_dist_and_k( Chris@16: function, k, p, k, &result, Policy()) Chris@16: && Chris@16: binomial_detail::check_dist_and_prob( Chris@16: function, k, p, alpha, &result, Policy())) Chris@16: { return result; } Chris@16: Chris@16: result = ibetac_invb(k + 1, p, alpha, Policy()); // returns n - k Chris@16: return result + k; Chris@16: } Chris@16: Chris@16: static RealType find_maximum_number_of_trials( Chris@16: RealType k, // number of events Chris@16: RealType p, // success fraction Chris@16: RealType alpha) // risk level Chris@16: { Chris@16: static const char* function = "boost::math::binomial_distribution<%1%>::find_maximum_number_of_trials"; Chris@16: // Error checks: Chris@16: RealType result = 0; Chris@16: if(false == binomial_detail::check_dist_and_k( Chris@16: function, k, p, k, &result, Policy()) Chris@16: && Chris@16: binomial_detail::check_dist_and_prob( Chris@16: function, k, p, alpha, &result, Policy())) Chris@16: { return result; } Chris@16: Chris@16: result = ibeta_invb(k + 1, p, alpha, Policy()); // returns n - k Chris@16: return result + k; Chris@16: } Chris@16: Chris@16: private: Chris@16: RealType m_n; // Not sure if this shouldn't be an int? Chris@16: RealType m_p; // success_fraction Chris@16: }; // template class binomial_distribution Chris@16: Chris@16: typedef binomial_distribution<> binomial; Chris@16: // typedef binomial_distribution binomial; Chris@16: // IS now included since no longer a name clash with function binomial. Chris@16: //typedef binomial_distribution binomial; // Reserved name of type double. Chris@16: Chris@16: template Chris@16: const std::pair range(const binomial_distribution& dist) Chris@16: { // Range of permissible values for random variable k. Chris@16: using boost::math::tools::max_value; Chris@16: return std::pair(static_cast(0), dist.trials()); Chris@16: } Chris@16: Chris@16: template Chris@16: const std::pair support(const binomial_distribution& 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: return std::pair(static_cast(0), dist.trials()); Chris@16: } Chris@16: Chris@16: template Chris@16: inline RealType mean(const binomial_distribution& dist) Chris@16: { // Mean of Binomial distribution = np. Chris@16: return dist.trials() * dist.success_fraction(); Chris@16: } // mean Chris@16: Chris@16: template Chris@16: inline RealType variance(const binomial_distribution& dist) Chris@16: { // Variance of Binomial distribution = np(1-p). Chris@16: return dist.trials() * dist.success_fraction() * (1 - dist.success_fraction()); Chris@16: } // variance Chris@16: Chris@16: template Chris@16: RealType pdf(const binomial_distribution& dist, const RealType& k) Chris@16: { // Probability Density/Mass Function. Chris@16: BOOST_FPU_EXCEPTION_GUARD Chris@16: Chris@16: BOOST_MATH_STD_USING // for ADL of std functions Chris@16: Chris@16: RealType n = dist.trials(); Chris@16: Chris@16: // Error check: Chris@16: RealType result = 0; // initialization silences some compiler warnings Chris@16: if(false == binomial_detail::check_dist_and_k( Chris@16: "boost::math::pdf(binomial_distribution<%1%> const&, %1%)", Chris@16: n, Chris@16: dist.success_fraction(), Chris@16: k, Chris@16: &result, Policy())) Chris@16: { Chris@16: return result; Chris@16: } Chris@16: Chris@16: // Special cases of success_fraction, regardless of k successes and regardless of n trials. Chris@16: if (dist.success_fraction() == 0) Chris@16: { // probability of zero successes is 1: Chris@16: return static_cast(k == 0 ? 1 : 0); Chris@16: } Chris@16: if (dist.success_fraction() == 1) Chris@16: { // probability of n successes is 1: Chris@16: return static_cast(k == n ? 1 : 0); Chris@16: } 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: if (n == 0) Chris@16: { Chris@16: return 1; // Probability = 1 = certainty. Chris@16: } Chris@16: if (k == 0) Chris@16: { // binomial coeffic (n 0) = 1, Chris@16: // n ^ 0 = 1 Chris@16: return pow(1 - dist.success_fraction(), n); Chris@16: } Chris@16: if (k == n) Chris@16: { // binomial coeffic (n n) = 1, Chris@16: // n ^ 0 = 1 Chris@16: return pow(dist.success_fraction(), k); // * pow((1 - dist.success_fraction()), (n - k)) = 1 Chris@16: } Chris@16: Chris@16: // Probability of getting exactly k successes Chris@16: // if C(n, k) is the binomial coefficient then: Chris@16: // Chris@16: // f(k; n,p) = C(n, k) * p^k * (1-p)^(n-k) Chris@16: // = (n!/(k!(n-k)!)) * p^k * (1-p)^(n-k) Chris@16: // = (tgamma(n+1) / (tgamma(k+1)*tgamma(n-k+1))) * p^k * (1-p)^(n-k) Chris@16: // = p^k (1-p)^(n-k) / (beta(k+1, n-k+1) * (n+1)) Chris@16: // = ibeta_derivative(k+1, n-k+1, p) / (n+1) Chris@16: // Chris@16: using boost::math::ibeta_derivative; // a, b, x Chris@16: return ibeta_derivative(k+1, n-k+1, dist.success_fraction(), Policy()) / (n+1); Chris@16: Chris@16: } // pdf Chris@16: Chris@16: template Chris@16: inline RealType cdf(const binomial_distribution& dist, const RealType& k) Chris@16: { // Cumulative Distribution Function Binomial. Chris@16: // The random variate k is the number of successes in n trials. 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: // Returns the sum of the terms 0 through k of the Binomial Probability Density/Mass: Chris@16: // Chris@16: // i=k Chris@16: // -- ( n ) i n-i Chris@16: // > | | p (1-p) Chris@16: // -- ( i ) Chris@16: // i=0 Chris@16: Chris@16: // The terms are not summed directly instead Chris@16: // the incomplete beta integral is employed, Chris@16: // according to the formula: Chris@16: // P = I[1-p]( n-k, k+1). Chris@16: // = 1 - I[p](k + 1, n - k) Chris@16: Chris@16: BOOST_MATH_STD_USING // for ADL of std functions Chris@16: Chris@16: RealType n = dist.trials(); Chris@16: RealType p = dist.success_fraction(); Chris@16: Chris@16: // Error check: Chris@16: RealType result = 0; Chris@16: if(false == binomial_detail::check_dist_and_k( Chris@16: "boost::math::cdf(binomial_distribution<%1%> const&, %1%)", Chris@16: n, Chris@16: p, Chris@16: k, Chris@16: &result, Policy())) Chris@16: { Chris@16: return result; Chris@16: } Chris@16: if (k == n) Chris@16: { Chris@16: return 1; Chris@16: } Chris@16: Chris@16: // Special cases, regardless of k. Chris@16: if (p == 0) Chris@16: { // This need explanation: Chris@16: // the pdf is zero for all cases except when k == 0. Chris@16: // For zero p the probability of zero successes is one. Chris@16: // Therefore the cdf is always 1: Chris@16: // the probability of k or *fewer* successes is always 1 Chris@16: // if there are never any successes! Chris@16: return 1; Chris@16: } Chris@16: if (p == 1) Chris@16: { // This is correct but needs explanation: Chris@16: // when k = 1 Chris@16: // all the cdf and pdf values are zero *except* when k == n, Chris@16: // and that case has been handled above already. Chris@16: return 0; Chris@16: } Chris@16: // Chris@16: // P = I[1-p](n - k, k + 1) Chris@16: // = 1 - I[p](k + 1, n - k) Chris@16: // Use of ibetac here prevents cancellation errors in calculating Chris@16: // 1-p if p is very small, perhaps smaller than machine epsilon. Chris@16: // Chris@16: // Note that we do not use a finite sum here, since the incomplete Chris@16: // beta uses a finite sum internally for integer arguments, so Chris@16: // we'll just let it take care of the necessary logic. Chris@16: // Chris@16: return ibetac(k + 1, n - k, p, Policy()); Chris@16: } // binomial cdf Chris@16: Chris@16: template Chris@16: inline RealType cdf(const complemented2_type, RealType>& c) Chris@16: { // Complemented Cumulative Distribution Function Binomial. Chris@16: // The random variate k is the number of successes in n trials. 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: // Returns the sum of the terms k+1 through n of the Binomial Probability Density/Mass: Chris@16: // Chris@16: // i=n Chris@16: // -- ( n ) i n-i Chris@16: // > | | p (1-p) Chris@16: // -- ( i ) Chris@16: // i=k+1 Chris@16: Chris@16: // The terms are not summed directly instead Chris@16: // the incomplete beta integral is employed, Chris@16: // according to the formula: Chris@16: // Q = 1 -I[1-p]( n-k, k+1). Chris@16: // = I[p](k + 1, n - k) Chris@16: Chris@16: BOOST_MATH_STD_USING // for ADL of std functions Chris@16: Chris@16: RealType const& k = c.param; Chris@16: binomial_distribution const& dist = c.dist; Chris@16: RealType n = dist.trials(); Chris@16: RealType p = dist.success_fraction(); Chris@16: Chris@16: // Error checks: Chris@16: RealType result = 0; Chris@16: if(false == binomial_detail::check_dist_and_k( Chris@16: "boost::math::cdf(binomial_distribution<%1%> const&, %1%)", Chris@16: n, Chris@16: p, Chris@16: k, Chris@16: &result, Policy())) Chris@16: { Chris@16: return result; Chris@16: } Chris@16: Chris@16: if (k == n) Chris@16: { // Probability of greater than n successes is necessarily zero: Chris@16: return 0; Chris@16: } Chris@16: Chris@16: // Special cases, regardless of k. Chris@16: if (p == 0) Chris@16: { Chris@16: // This need explanation: the pdf is zero for all Chris@16: // cases except when k == 0. For zero p the probability Chris@16: // of zero successes is one. Therefore the cdf is always Chris@16: // 1: the probability of *more than* k successes is always 0 Chris@16: // if there are never any successes! Chris@16: return 0; Chris@16: } Chris@16: if (p == 1) Chris@16: { Chris@16: // This needs explanation, when p = 1 Chris@16: // we always have n successes, so the probability Chris@16: // of more than k successes is 1 as long as k < n. Chris@16: // The k == n case has already been handled above. Chris@16: return 1; Chris@16: } Chris@16: // Chris@16: // Calculate cdf binomial using the incomplete beta function. Chris@16: // Q = 1 -I[1-p](n - k, k + 1) Chris@16: // = I[p](k + 1, n - k) Chris@16: // Use of ibeta here prevents cancellation errors in calculating Chris@16: // 1-p if p is very small, perhaps smaller than machine epsilon. Chris@16: // Chris@16: // Note that we do not use a finite sum here, since the incomplete Chris@16: // beta uses a finite sum internally for integer arguments, so Chris@16: // we'll just let it take care of the necessary logic. Chris@16: // Chris@16: return ibeta(k + 1, n - k, p, Policy()); Chris@16: } // binomial cdf Chris@16: Chris@16: template Chris@16: inline RealType quantile(const binomial_distribution& dist, const RealType& p) Chris@16: { Chris@16: return binomial_detail::quantile_imp(dist, p, RealType(1-p), false); Chris@16: } // quantile Chris@16: Chris@16: template Chris@16: RealType quantile(const complemented2_type, RealType>& c) Chris@16: { Chris@16: return binomial_detail::quantile_imp(c.dist, RealType(1-c.param), c.param, true); Chris@16: } // quantile Chris@16: Chris@16: template Chris@16: inline RealType mode(const binomial_distribution& dist) Chris@16: { Chris@16: BOOST_MATH_STD_USING // ADL of std functions. Chris@16: RealType p = dist.success_fraction(); Chris@16: RealType n = dist.trials(); Chris@16: return floor(p * (n + 1)); Chris@16: } Chris@16: Chris@16: template Chris@16: inline RealType median(const binomial_distribution& dist) Chris@16: { // Bounds for the median of the negative binomial distribution Chris@16: // VAN DE VEN R. ; WEBER N. C. ; Chris@16: // Univ. Sydney, school mathematics statistics, Sydney N.S.W. 2006, AUSTRALIE Chris@16: // Metrika (Metrika) ISSN 0026-1335 CODEN MTRKA8 Chris@16: // 1993, vol. 40, no3-4, pp. 185-189 (4 ref.) Chris@16: Chris@16: // Bounds for median and 50 percetage point of binomial and negative binomial distribution Chris@16: // Metrika, ISSN 0026-1335 (Print) 1435-926X (Online) Chris@16: // Volume 41, Number 1 / December, 1994, DOI 10.1007/BF01895303 Chris@16: BOOST_MATH_STD_USING // ADL of std functions. Chris@16: RealType p = dist.success_fraction(); Chris@16: RealType n = dist.trials(); Chris@16: // Wikipedia says one of floor(np) -1, floor (np), floor(np) +1 Chris@16: return floor(p * n); // Chose the middle value. Chris@16: } Chris@16: Chris@16: template Chris@16: inline RealType skewness(const binomial_distribution& dist) Chris@16: { Chris@16: BOOST_MATH_STD_USING // ADL of std functions. Chris@16: RealType p = dist.success_fraction(); Chris@16: RealType n = dist.trials(); Chris@16: return (1 - 2 * p) / sqrt(n * p * (1 - p)); Chris@16: } Chris@16: Chris@16: template Chris@16: inline RealType kurtosis(const binomial_distribution& dist) Chris@16: { Chris@16: RealType p = dist.success_fraction(); Chris@16: RealType n = dist.trials(); Chris@16: return 3 - 6 / n + 1 / (n * p * (1 - p)); Chris@16: } Chris@16: Chris@16: template Chris@16: inline RealType kurtosis_excess(const binomial_distribution& dist) Chris@16: { Chris@16: RealType p = dist.success_fraction(); Chris@16: RealType q = 1 - p; Chris@16: RealType n = dist.trials(); Chris@16: return (1 - 6 * p * q) / (n * p * q); Chris@16: } 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 Chris@16: Chris@16: #endif // BOOST_MATH_SPECIAL_BINOMIAL_HPP Chris@16: Chris@16: