Chris@16: // boost\math\distributions\beta.hpp Chris@16: Chris@16: // Copyright John Maddock 2006. Chris@16: // Copyright Paul A. Bristow 2006. 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/Beta_distribution Chris@16: // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366h.htm Chris@16: // http://mathworld.wolfram.com/BetaDistribution.html Chris@16: Chris@16: // The Beta Distribution is a continuous probability distribution. Chris@16: // The beta distribution is used to model events which are constrained to take place Chris@16: // within an interval defined by maxima and minima, Chris@16: // so is used extensively in PERT and other project management systems Chris@16: // to describe the time to completion. Chris@16: // The cdf of the beta distribution is used as a convenient way Chris@16: // of obtaining the sum over a set of binomial outcomes. Chris@16: // The beta distribution is also used in Bayesian statistics. Chris@16: Chris@16: #ifndef BOOST_MATH_DIST_BETA_HPP Chris@16: #define BOOST_MATH_DIST_BETA_HPP Chris@16: Chris@16: #include Chris@16: #include // for beta. Chris@16: #include // complements. Chris@16: #include // error checks Chris@16: #include // isnan. Chris@16: #include // for root finding. Chris@16: Chris@16: #if defined (BOOST_MSVC) Chris@16: # pragma warning(push) Chris@16: # pragma warning(disable: 4702) // unreachable code Chris@16: // in domain_error_imp in error_handling Chris@16: #endif Chris@16: Chris@16: #include Chris@16: Chris@16: namespace boost Chris@16: { Chris@16: namespace math Chris@16: { Chris@16: namespace beta_detail Chris@16: { Chris@16: // Common error checking routines for beta distribution functions: Chris@16: template Chris@16: inline bool check_alpha(const char* function, const RealType& alpha, RealType* result, const Policy& pol) Chris@16: { Chris@16: if(!(boost::math::isfinite)(alpha) || (alpha <= 0)) Chris@16: { Chris@16: *result = policies::raise_domain_error( Chris@16: function, Chris@16: "Alpha argument is %1%, but must be > 0 !", alpha, pol); Chris@16: return false; Chris@16: } Chris@16: return true; Chris@16: } // bool check_alpha Chris@16: Chris@16: template Chris@16: inline bool check_beta(const char* function, const RealType& beta, RealType* result, const Policy& pol) Chris@16: { Chris@16: if(!(boost::math::isfinite)(beta) || (beta <= 0)) Chris@16: { Chris@16: *result = policies::raise_domain_error( Chris@16: function, Chris@16: "Beta argument is %1%, but must be > 0 !", beta, pol); Chris@16: return false; Chris@16: } Chris@16: return true; Chris@16: } // bool check_beta Chris@16: Chris@16: template Chris@16: inline bool check_prob(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: "Probability argument is %1%, but must be >= 0 and <= 1 !", p, pol); Chris@16: return false; Chris@16: } Chris@16: return true; Chris@16: } // bool check_prob Chris@16: Chris@16: template Chris@16: inline bool check_x(const char* function, const RealType& x, RealType* result, const Policy& pol) Chris@16: { Chris@16: if(!(boost::math::isfinite)(x) || (x < 0) || (x > 1)) Chris@16: { Chris@16: *result = policies::raise_domain_error( Chris@16: function, Chris@16: "x argument is %1%, but must be >= 0 and <= 1 !", x, pol); Chris@16: return false; Chris@16: } Chris@16: return true; Chris@16: } // bool check_x Chris@16: Chris@16: template Chris@16: inline bool check_dist(const char* function, const RealType& alpha, const RealType& beta, RealType* result, const Policy& pol) Chris@16: { // Check both alpha and beta. Chris@16: return check_alpha(function, alpha, result, pol) Chris@16: && check_beta(function, beta, result, pol); Chris@16: } // bool check_dist Chris@16: Chris@16: template Chris@16: inline bool check_dist_and_x(const char* function, const RealType& alpha, const RealType& beta, RealType x, RealType* result, const Policy& pol) Chris@16: { Chris@16: return check_dist(function, alpha, beta, result, pol) Chris@16: && beta_detail::check_x(function, x, result, pol); Chris@16: } // bool check_dist_and_x Chris@16: Chris@16: template Chris@16: inline bool check_dist_and_prob(const char* function, const RealType& alpha, const RealType& beta, RealType p, RealType* result, const Policy& pol) Chris@16: { Chris@16: return check_dist(function, alpha, beta, result, pol) Chris@16: && check_prob(function, p, result, pol); Chris@16: } // bool check_dist_and_prob Chris@16: Chris@16: template Chris@16: inline bool check_mean(const char* function, const RealType& mean, RealType* result, const Policy& pol) Chris@16: { Chris@16: if(!(boost::math::isfinite)(mean) || (mean <= 0)) Chris@16: { Chris@16: *result = policies::raise_domain_error( Chris@16: function, Chris@16: "mean argument is %1%, but must be > 0 !", mean, pol); Chris@16: return false; Chris@16: } Chris@16: return true; Chris@16: } // bool check_mean Chris@16: template Chris@16: inline bool check_variance(const char* function, const RealType& variance, RealType* result, const Policy& pol) Chris@16: { Chris@16: if(!(boost::math::isfinite)(variance) || (variance <= 0)) Chris@16: { Chris@16: *result = policies::raise_domain_error( Chris@16: function, Chris@16: "variance argument is %1%, but must be > 0 !", variance, pol); Chris@16: return false; Chris@16: } Chris@16: return true; Chris@16: } // bool check_variance Chris@16: } // namespace beta_detail Chris@16: Chris@16: // typedef beta_distribution beta; Chris@16: // is deliberately NOT included to avoid a name clash with the beta function. Chris@16: // Use beta_distribution<> mybeta(...) to construct type double. Chris@16: Chris@16: template > Chris@16: class beta_distribution Chris@16: { Chris@16: public: Chris@16: typedef RealType value_type; Chris@16: typedef Policy policy_type; Chris@16: Chris@16: beta_distribution(RealType l_alpha = 1, RealType l_beta = 1) : m_alpha(l_alpha), m_beta(l_beta) Chris@16: { Chris@16: RealType result; Chris@16: beta_detail::check_dist( Chris@16: "boost::math::beta_distribution<%1%>::beta_distribution", Chris@16: m_alpha, Chris@16: m_beta, Chris@16: &result, Policy()); Chris@16: } // beta_distribution constructor. Chris@16: // Accessor functions: Chris@16: RealType alpha() const Chris@16: { Chris@16: return m_alpha; Chris@16: } Chris@16: RealType beta() const Chris@16: { // . Chris@16: return m_beta; Chris@16: } Chris@16: Chris@16: // Estimation of the alpha & beta parameters. Chris@16: // http://en.wikipedia.org/wiki/Beta_distribution Chris@16: // gives formulae in section on parameter estimation. Chris@16: // Also NIST EDA page 3 & 4 give the same. Chris@16: // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366h.htm Chris@16: // http://www.epi.ucdavis.edu/diagnostictests/betabuster.html Chris@16: Chris@16: static RealType find_alpha( Chris@16: RealType mean, // Expected value of mean. Chris@16: RealType variance) // Expected value of variance. Chris@16: { Chris@16: static const char* function = "boost::math::beta_distribution<%1%>::find_alpha"; Chris@16: RealType result = 0; // of error checks. Chris@16: if(false == Chris@16: ( Chris@16: beta_detail::check_mean(function, mean, &result, Policy()) Chris@16: && beta_detail::check_variance(function, variance, &result, Policy()) Chris@16: ) Chris@16: ) Chris@16: { Chris@16: return result; Chris@16: } Chris@16: return mean * (( (mean * (1 - mean)) / variance)- 1); Chris@16: } // RealType find_alpha Chris@16: Chris@16: static RealType find_beta( Chris@16: RealType mean, // Expected value of mean. Chris@16: RealType variance) // Expected value of variance. Chris@16: { Chris@16: static const char* function = "boost::math::beta_distribution<%1%>::find_beta"; Chris@16: RealType result = 0; // of error checks. Chris@16: if(false == Chris@16: ( Chris@16: beta_detail::check_mean(function, mean, &result, Policy()) Chris@16: && Chris@16: beta_detail::check_variance(function, variance, &result, Policy()) Chris@16: ) Chris@16: ) Chris@16: { Chris@16: return result; Chris@16: } Chris@16: return (1 - mean) * (((mean * (1 - mean)) /variance)-1); Chris@16: } // RealType find_beta Chris@16: Chris@16: // Estimate alpha & beta from either alpha or beta, and x and probability. Chris@16: // Uses for these parameter estimators are unclear. Chris@16: Chris@16: static RealType find_alpha( Chris@16: RealType beta, // from beta. Chris@16: RealType x, // x. Chris@16: RealType probability) // cdf Chris@16: { Chris@16: static const char* function = "boost::math::beta_distribution<%1%>::find_alpha"; Chris@16: RealType result = 0; // of error checks. Chris@16: if(false == Chris@16: ( Chris@16: beta_detail::check_prob(function, probability, &result, Policy()) Chris@16: && Chris@16: beta_detail::check_beta(function, beta, &result, Policy()) Chris@16: && Chris@16: beta_detail::check_x(function, x, &result, Policy()) Chris@16: ) Chris@16: ) Chris@16: { Chris@16: return result; Chris@16: } Chris@16: return ibeta_inva(beta, x, probability, Policy()); Chris@16: } // RealType find_alpha(beta, a, probability) Chris@16: Chris@16: static RealType find_beta( Chris@16: // ibeta_invb(T b, T x, T p); (alpha, x, cdf,) Chris@16: RealType alpha, // alpha. Chris@16: RealType x, // probability x. Chris@16: RealType probability) // probability cdf. Chris@16: { Chris@16: static const char* function = "boost::math::beta_distribution<%1%>::find_beta"; Chris@16: RealType result = 0; // of error checks. Chris@16: if(false == Chris@16: ( Chris@16: beta_detail::check_prob(function, probability, &result, Policy()) Chris@16: && Chris@16: beta_detail::check_alpha(function, alpha, &result, Policy()) Chris@16: && Chris@16: beta_detail::check_x(function, x, &result, Policy()) Chris@16: ) Chris@16: ) Chris@16: { Chris@16: return result; Chris@16: } Chris@16: return ibeta_invb(alpha, x, probability, Policy()); Chris@16: } // RealType find_beta(alpha, x, probability) Chris@16: Chris@16: private: Chris@16: RealType m_alpha; // Two parameters of the beta distribution. Chris@16: RealType m_beta; Chris@16: }; // template class beta_distribution Chris@16: Chris@16: template Chris@16: inline const std::pair range(const beta_distribution& /* dist */) Chris@16: { // Range of permissible values for random variable x. Chris@16: using boost::math::tools::max_value; Chris@16: return std::pair(static_cast(0), static_cast(1)); Chris@16: } Chris@16: Chris@16: template Chris@16: inline const std::pair support(const beta_distribution& /* dist */) Chris@16: { // Range of supported values for random variable x. 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), static_cast(1)); Chris@16: } Chris@16: Chris@16: template Chris@16: inline RealType mean(const beta_distribution& dist) Chris@16: { // Mean of beta distribution = np. Chris@16: return dist.alpha() / (dist.alpha() + dist.beta()); Chris@16: } // mean Chris@16: Chris@16: template Chris@16: inline RealType variance(const beta_distribution& dist) Chris@16: { // Variance of beta distribution = np(1-p). Chris@16: RealType a = dist.alpha(); Chris@16: RealType b = dist.beta(); Chris@16: return (a * b) / ((a + b ) * (a + b) * (a + b + 1)); Chris@16: } // variance Chris@16: Chris@16: template Chris@16: inline RealType mode(const beta_distribution& dist) Chris@16: { Chris@16: static const char* function = "boost::math::mode(beta_distribution<%1%> const&)"; Chris@16: Chris@16: RealType result; Chris@16: if ((dist.alpha() <= 1)) Chris@16: { Chris@16: result = policies::raise_domain_error( Chris@16: function, Chris@16: "mode undefined for alpha = %1%, must be > 1!", dist.alpha(), Policy()); Chris@16: return result; Chris@16: } Chris@16: Chris@16: if ((dist.beta() <= 1)) Chris@16: { Chris@16: result = policies::raise_domain_error( Chris@16: function, Chris@16: "mode undefined for beta = %1%, must be > 1!", dist.beta(), Policy()); Chris@16: return result; Chris@16: } Chris@16: RealType a = dist.alpha(); Chris@16: RealType b = dist.beta(); Chris@16: return (a-1) / (a + b - 2); Chris@16: } // mode Chris@16: Chris@16: //template Chris@16: //inline RealType median(const beta_distribution& dist) Chris@16: //{ // Median of beta distribution is not defined. Chris@16: // return tools::domain_error(function, "Median is not implemented, result is %1%!", std::numeric_limits::quiet_NaN()); Chris@16: //} // median Chris@16: Chris@16: //But WILL be provided by the derived accessor as quantile(0.5). Chris@16: Chris@16: template Chris@16: inline RealType skewness(const beta_distribution& dist) Chris@16: { Chris@16: BOOST_MATH_STD_USING // ADL of std functions. Chris@16: RealType a = dist.alpha(); Chris@16: RealType b = dist.beta(); Chris@16: return (2 * (b-a) * sqrt(a + b + 1)) / ((a + b + 2) * sqrt(a * b)); Chris@16: } // skewness Chris@16: Chris@16: template Chris@16: inline RealType kurtosis_excess(const beta_distribution& dist) Chris@16: { Chris@16: RealType a = dist.alpha(); Chris@16: RealType b = dist.beta(); Chris@16: RealType a_2 = a * a; Chris@16: RealType n = 6 * (a_2 * a - a_2 * (2 * b - 1) + b * b * (b + 1) - 2 * a * b * (b + 2)); Chris@16: RealType d = a * b * (a + b + 2) * (a + b + 3); Chris@16: return n / d; Chris@16: } // kurtosis_excess Chris@16: Chris@16: template Chris@16: inline RealType kurtosis(const beta_distribution& dist) Chris@16: { Chris@16: return 3 + kurtosis_excess(dist); Chris@16: } // kurtosis Chris@16: Chris@16: template Chris@16: inline RealType pdf(const beta_distribution& dist, const RealType& x) Chris@16: { // Probability Density/Mass Function. Chris@16: BOOST_FPU_EXCEPTION_GUARD Chris@16: Chris@16: static const char* function = "boost::math::pdf(beta_distribution<%1%> const&, %1%)"; Chris@16: Chris@16: BOOST_MATH_STD_USING // for ADL of std functions Chris@16: Chris@16: RealType a = dist.alpha(); Chris@16: RealType b = dist.beta(); Chris@16: Chris@16: // Argument checks: Chris@16: RealType result = 0; Chris@16: if(false == beta_detail::check_dist_and_x( Chris@16: function, Chris@16: a, b, x, Chris@16: &result, Policy())) Chris@16: { Chris@16: return result; Chris@16: } Chris@16: using boost::math::beta; Chris@16: return ibeta_derivative(a, b, x, Policy()); Chris@16: } // pdf Chris@16: Chris@16: template Chris@16: inline RealType cdf(const beta_distribution& dist, const RealType& x) Chris@16: { // Cumulative Distribution Function beta. Chris@16: BOOST_MATH_STD_USING // for ADL of std functions Chris@16: Chris@16: static const char* function = "boost::math::cdf(beta_distribution<%1%> const&, %1%)"; Chris@16: Chris@16: RealType a = dist.alpha(); Chris@16: RealType b = dist.beta(); Chris@16: Chris@16: // Argument checks: Chris@16: RealType result = 0; Chris@16: if(false == beta_detail::check_dist_and_x( Chris@16: function, Chris@16: a, b, x, Chris@16: &result, Policy())) Chris@16: { Chris@16: return result; Chris@16: } Chris@16: // Special cases: Chris@16: if (x == 0) Chris@16: { Chris@16: return 0; Chris@16: } Chris@16: else if (x == 1) Chris@16: { Chris@16: return 1; Chris@16: } Chris@16: return ibeta(a, b, x, Policy()); Chris@16: } // beta cdf Chris@16: Chris@16: template Chris@16: inline RealType cdf(const complemented2_type, RealType>& c) Chris@16: { // Complemented Cumulative Distribution Function beta. Chris@16: Chris@16: BOOST_MATH_STD_USING // for ADL of std functions Chris@16: Chris@16: static const char* function = "boost::math::cdf(beta_distribution<%1%> const&, %1%)"; Chris@16: Chris@16: RealType const& x = c.param; Chris@16: beta_distribution const& dist = c.dist; Chris@16: RealType a = dist.alpha(); Chris@16: RealType b = dist.beta(); Chris@16: Chris@16: // Argument checks: Chris@16: RealType result = 0; Chris@16: if(false == beta_detail::check_dist_and_x( Chris@16: function, Chris@16: a, b, x, Chris@16: &result, Policy())) Chris@16: { Chris@16: return result; Chris@16: } Chris@16: if (x == 0) Chris@16: { Chris@16: return 1; Chris@16: } Chris@16: else if (x == 1) Chris@16: { Chris@16: return 0; Chris@16: } Chris@16: // Calculate cdf beta using the incomplete beta function. Chris@16: // Use of ibeta here prevents cancellation errors in calculating Chris@16: // 1 - x if x is very small, perhaps smaller than machine epsilon. Chris@16: return ibetac(a, b, x, Policy()); Chris@16: } // beta cdf Chris@16: Chris@16: template Chris@16: inline RealType quantile(const beta_distribution& dist, const RealType& p) Chris@16: { // Quantile or Percent Point beta function or Chris@16: // Inverse Cumulative probability distribution function CDF. Chris@16: // Return x (0 <= x <= 1), Chris@16: // for a given probability p (0 <= p <= 1). Chris@16: // These functions take a probability as an argument Chris@16: // and return a value such that the probability that a random variable x Chris@16: // will be less than or equal to that value Chris@16: // is whatever probability you supplied as an argument. Chris@16: Chris@16: static const char* function = "boost::math::quantile(beta_distribution<%1%> const&, %1%)"; Chris@16: Chris@16: RealType result = 0; // of argument checks: Chris@16: RealType a = dist.alpha(); Chris@16: RealType b = dist.beta(); Chris@16: if(false == beta_detail::check_dist_and_prob( Chris@16: function, Chris@16: a, b, p, Chris@16: &result, Policy())) Chris@16: { Chris@16: return result; Chris@16: } Chris@16: // Special cases: Chris@16: if (p == 0) Chris@16: { Chris@16: return 0; Chris@16: } Chris@16: if (p == 1) Chris@16: { Chris@16: return 1; Chris@16: } Chris@16: return ibeta_inv(a, b, p, static_cast(0), Policy()); Chris@16: } // quantile Chris@16: Chris@16: template Chris@16: inline RealType quantile(const complemented2_type, RealType>& c) Chris@16: { // Complement Quantile or Percent Point beta function . Chris@16: // Return the number of expected x for a given Chris@16: // complement of the probability q. Chris@16: Chris@16: static const char* function = "boost::math::quantile(beta_distribution<%1%> const&, %1%)"; Chris@16: Chris@16: // Chris@16: // Error checks: Chris@16: RealType q = c.param; Chris@16: const beta_distribution& dist = c.dist; Chris@16: RealType result = 0; Chris@16: RealType a = dist.alpha(); Chris@16: RealType b = dist.beta(); Chris@16: if(false == beta_detail::check_dist_and_prob( Chris@16: function, Chris@16: a, Chris@16: b, Chris@16: q, Chris@16: &result, Policy())) Chris@16: { Chris@16: return result; Chris@16: } Chris@16: // Special cases: Chris@16: if(q == 1) Chris@16: { Chris@16: return 0; Chris@16: } Chris@16: if(q == 0) Chris@16: { Chris@16: return 1; Chris@16: } Chris@16: Chris@16: return ibetac_inv(a, b, q, static_cast(0), Policy()); Chris@16: } // Quantile Complement Chris@16: Chris@16: } // namespace math Chris@16: } // namespace boost Chris@16: Chris@16: // This include must be at the end, *after* the accessors Chris@16: // for this distribution have been defined, in order to Chris@16: // keep compilers that support two-phase lookup happy. Chris@16: #include Chris@16: Chris@16: #if defined (BOOST_MSVC) Chris@16: # pragma warning(pop) Chris@16: #endif Chris@16: Chris@16: #endif // BOOST_MATH_DIST_BETA_HPP Chris@16: Chris@16: