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