585 lines · plain
1// boost\math\distributions\beta.hpp2 3// Copyright John Maddock 2006.4// Copyright Paul A. Bristow 2006.5// Copyright Matt Borland 2023.6 7// Use, modification and distribution are subject to the8// Boost Software License, Version 1.0.9// (See accompanying file LICENSE_1_0.txt10// or copy at http://www.boost.org/LICENSE_1_0.txt)11 12// http://en.wikipedia.org/wiki/Beta_distribution13// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366h.htm14// http://mathworld.wolfram.com/BetaDistribution.html15 16// The Beta Distribution is a continuous probability distribution.17// The beta distribution is used to model events which are constrained to take place18// within an interval defined by maxima and minima,19// so is used extensively in PERT and other project management systems20// to describe the time to completion.21// The cdf of the beta distribution is used as a convenient way22// of obtaining the sum over a set of binomial outcomes.23// The beta distribution is also used in Bayesian statistics.24 25#ifndef BOOST_MATH_DIST_BETA_HPP26#define BOOST_MATH_DIST_BETA_HPP27 28#include <boost/math/tools/config.hpp>29#include <boost/math/tools/tuple.hpp>30#include <boost/math/distributions/fwd.hpp>31#include <boost/math/special_functions/beta.hpp> // for beta.32#include <boost/math/distributions/complement.hpp> // complements.33#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks34#include <boost/math/special_functions/fpclassify.hpp> // isnan.35#include <boost/math/tools/roots.hpp> // for root finding.36#include <boost/math/policies/error_handling.hpp>37 38#if defined (BOOST_MSVC)39# pragma warning(push)40# pragma warning(disable: 4702) // unreachable code41// in domain_error_imp in error_handling42#endif43 44namespace boost45{46 namespace math47 {48 namespace beta_detail49 {50 // Common error checking routines for beta distribution functions:51 template <class RealType, class Policy>52 BOOST_MATH_GPU_ENABLED inline bool check_alpha(const char* function, const RealType& alpha, RealType* result, const Policy& pol)53 {54 if(!(boost::math::isfinite)(alpha) || (alpha <= 0))55 {56 *result = policies::raise_domain_error<RealType>(57 function,58 "Alpha argument is %1%, but must be > 0 !", alpha, pol);59 return false;60 }61 return true;62 } // bool check_alpha63 64 template <class RealType, class Policy>65 BOOST_MATH_GPU_ENABLED inline bool check_beta(const char* function, const RealType& beta, RealType* result, const Policy& pol)66 {67 if(!(boost::math::isfinite)(beta) || (beta <= 0))68 {69 *result = policies::raise_domain_error<RealType>(70 function,71 "Beta argument is %1%, but must be > 0 !", beta, pol);72 return false;73 }74 return true;75 } // bool check_beta76 77 template <class RealType, class Policy>78 BOOST_MATH_GPU_ENABLED inline bool check_prob(const char* function, const RealType& p, RealType* result, const Policy& pol)79 {80 if((p < 0) || (p > 1) || !(boost::math::isfinite)(p))81 {82 *result = policies::raise_domain_error<RealType>(83 function,84 "Probability argument is %1%, but must be >= 0 and <= 1 !", p, pol);85 return false;86 }87 return true;88 } // bool check_prob89 90 template <class RealType, class Policy>91 BOOST_MATH_GPU_ENABLED inline bool check_x(const char* function, const RealType& x, RealType* result, const Policy& pol)92 {93 if(!(boost::math::isfinite)(x) || (x < 0) || (x > 1))94 {95 *result = policies::raise_domain_error<RealType>(96 function,97 "x argument is %1%, but must be >= 0 and <= 1 !", x, pol);98 return false;99 }100 return true;101 } // bool check_x102 103 template <class RealType, class Policy>104 BOOST_MATH_GPU_ENABLED inline bool check_dist(const char* function, const RealType& alpha, const RealType& beta, RealType* result, const Policy& pol)105 { // Check both alpha and beta.106 return check_alpha(function, alpha, result, pol)107 && check_beta(function, beta, result, pol);108 } // bool check_dist109 110 template <class RealType, class Policy>111 BOOST_MATH_GPU_ENABLED inline bool check_dist_and_x(const char* function, const RealType& alpha, const RealType& beta, RealType x, RealType* result, const Policy& pol)112 {113 return check_dist(function, alpha, beta, result, pol)114 && beta_detail::check_x(function, x, result, pol);115 } // bool check_dist_and_x116 117 template <class RealType, class Policy>118 BOOST_MATH_GPU_ENABLED inline bool check_dist_and_prob(const char* function, const RealType& alpha, const RealType& beta, RealType p, RealType* result, const Policy& pol)119 {120 return check_dist(function, alpha, beta, result, pol)121 && check_prob(function, p, result, pol);122 } // bool check_dist_and_prob123 124 template <class RealType, class Policy>125 BOOST_MATH_GPU_ENABLED inline bool check_mean(const char* function, const RealType& mean, RealType* result, const Policy& pol)126 {127 if(!(boost::math::isfinite)(mean) || (mean <= 0))128 {129 *result = policies::raise_domain_error<RealType>(130 function,131 "mean argument is %1%, but must be > 0 !", mean, pol);132 return false;133 }134 return true;135 } // bool check_mean136 template <class RealType, class Policy>137 BOOST_MATH_GPU_ENABLED inline bool check_variance(const char* function, const RealType& variance, RealType* result, const Policy& pol)138 {139 if(!(boost::math::isfinite)(variance) || (variance <= 0))140 {141 *result = policies::raise_domain_error<RealType>(142 function,143 "variance argument is %1%, but must be > 0 !", variance, pol);144 return false;145 }146 return true;147 } // bool check_variance148 } // namespace beta_detail149 150 // typedef beta_distribution<double> beta;151 // is deliberately NOT included to avoid a name clash with the beta function.152 // Use beta_distribution<> mybeta(...) to construct type double.153 154 template <class RealType = double, class Policy = policies::policy<> >155 class beta_distribution156 {157 public:158 typedef RealType value_type;159 typedef Policy policy_type;160 161 BOOST_MATH_GPU_ENABLED beta_distribution(RealType l_alpha = 1, RealType l_beta = 1) : m_alpha(l_alpha), m_beta(l_beta)162 {163 RealType result;164 beta_detail::check_dist(165 "boost::math::beta_distribution<%1%>::beta_distribution",166 m_alpha,167 m_beta,168 &result, Policy());169 } // beta_distribution constructor.170 // Accessor functions:171 BOOST_MATH_GPU_ENABLED RealType alpha() const172 {173 return m_alpha;174 }175 BOOST_MATH_GPU_ENABLED RealType beta() const176 { // .177 return m_beta;178 }179 180 // Estimation of the alpha & beta parameters.181 // http://en.wikipedia.org/wiki/Beta_distribution182 // gives formulae in section on parameter estimation.183 // Also NIST EDA page 3 & 4 give the same.184 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366h.htm185 // http://www.epi.ucdavis.edu/diagnostictests/betabuster.html186 187 BOOST_MATH_GPU_ENABLED static RealType find_alpha(188 RealType mean, // Expected value of mean.189 RealType variance) // Expected value of variance.190 {191 constexpr auto function = "boost::math::beta_distribution<%1%>::find_alpha";192 RealType result = 0; // of error checks.193 if(false ==194 (195 beta_detail::check_mean(function, mean, &result, Policy())196 && beta_detail::check_variance(function, variance, &result, Policy())197 )198 )199 {200 return result;201 }202 return mean * (( (mean * (1 - mean)) / variance)- 1);203 } // RealType find_alpha204 205 BOOST_MATH_GPU_ENABLED static RealType find_beta(206 RealType mean, // Expected value of mean.207 RealType variance) // Expected value of variance.208 {209 constexpr auto function = "boost::math::beta_distribution<%1%>::find_beta";210 RealType result = 0; // of error checks.211 if(false ==212 (213 beta_detail::check_mean(function, mean, &result, Policy())214 &&215 beta_detail::check_variance(function, variance, &result, Policy())216 )217 )218 {219 return result;220 }221 return (1 - mean) * (((mean * (1 - mean)) /variance)-1);222 } // RealType find_beta223 224 // Estimate alpha & beta from either alpha or beta, and x and probability.225 // Uses for these parameter estimators are unclear.226 227 BOOST_MATH_GPU_ENABLED static RealType find_alpha(228 RealType beta, // from beta.229 RealType x, // x.230 RealType probability) // cdf231 {232 constexpr auto function = "boost::math::beta_distribution<%1%>::find_alpha";233 RealType result = 0; // of error checks.234 if(false ==235 (236 beta_detail::check_prob(function, probability, &result, Policy())237 &&238 beta_detail::check_beta(function, beta, &result, Policy())239 &&240 beta_detail::check_x(function, x, &result, Policy())241 )242 )243 {244 return result;245 }246 return static_cast<RealType>(ibeta_inva(beta, x, probability, Policy()));247 } // RealType find_alpha(beta, a, probability)248 249 BOOST_MATH_GPU_ENABLED static RealType find_beta(250 // ibeta_invb(T b, T x, T p); (alpha, x, cdf,)251 RealType alpha, // alpha.252 RealType x, // probability x.253 RealType probability) // probability cdf.254 {255 constexpr auto function = "boost::math::beta_distribution<%1%>::find_beta";256 RealType result = 0; // of error checks.257 if(false ==258 (259 beta_detail::check_prob(function, probability, &result, Policy())260 &&261 beta_detail::check_alpha(function, alpha, &result, Policy())262 &&263 beta_detail::check_x(function, x, &result, Policy())264 )265 )266 {267 return result;268 }269 return static_cast<RealType>(ibeta_invb(alpha, x, probability, Policy()));270 } // RealType find_beta(alpha, x, probability)271 272 private:273 RealType m_alpha; // Two parameters of the beta distribution.274 RealType m_beta;275 }; // template <class RealType, class Policy> class beta_distribution276 277 #ifdef __cpp_deduction_guides278 template <class RealType>279 beta_distribution(RealType)->beta_distribution<typename boost::math::tools::promote_args<RealType>::type>;280 template <class RealType>281 beta_distribution(RealType, RealType)->beta_distribution<typename boost::math::tools::promote_args<RealType>::type>;282 #endif283 284 template <class RealType, class Policy>285 BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> range(const beta_distribution<RealType, Policy>& /* dist */)286 { // Range of permissible values for random variable x.287 using boost::math::tools::max_value;288 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1));289 }290 291 template <class RealType, class Policy>292 BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> support(const beta_distribution<RealType, Policy>& /* dist */)293 { // Range of supported values for random variable x.294 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.295 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1));296 }297 298 template <class RealType, class Policy>299 BOOST_MATH_GPU_ENABLED inline RealType mean(const beta_distribution<RealType, Policy>& dist)300 { // Mean of beta distribution = np.301 return dist.alpha() / (dist.alpha() + dist.beta());302 } // mean303 304 template <class RealType, class Policy>305 BOOST_MATH_GPU_ENABLED inline RealType variance(const beta_distribution<RealType, Policy>& dist)306 { // Variance of beta distribution = np(1-p).307 RealType a = dist.alpha();308 RealType b = dist.beta();309 return (a * b) / ((a + b ) * (a + b) * (a + b + 1));310 } // variance311 312 template <class RealType, class Policy>313 BOOST_MATH_GPU_ENABLED inline RealType mode(const beta_distribution<RealType, Policy>& dist)314 {315 constexpr auto function = "boost::math::mode(beta_distribution<%1%> const&)";316 317 RealType result;318 if ((dist.alpha() <= 1))319 {320 result = policies::raise_domain_error<RealType>(321 function,322 "mode undefined for alpha = %1%, must be > 1!", dist.alpha(), Policy());323 return result;324 }325 326 if ((dist.beta() <= 1))327 {328 result = policies::raise_domain_error<RealType>(329 function,330 "mode undefined for beta = %1%, must be > 1!", dist.beta(), Policy());331 return result;332 }333 RealType a = dist.alpha();334 RealType b = dist.beta();335 return (a-1) / (a + b - 2);336 } // mode337 338 //template <class RealType, class Policy>339 //inline RealType median(const beta_distribution<RealType, Policy>& dist)340 //{ // Median of beta distribution is not defined.341 // return tools::domain_error<RealType>(function, "Median is not implemented, result is %1%!", std::numeric_limits<RealType>::quiet_NaN());342 //} // median343 344 //But WILL be provided by the derived accessor as quantile(0.5).345 346 template <class RealType, class Policy>347 BOOST_MATH_GPU_ENABLED inline RealType skewness(const beta_distribution<RealType, Policy>& dist)348 {349 BOOST_MATH_STD_USING // ADL of std functions.350 RealType a = dist.alpha();351 RealType b = dist.beta();352 return (2 * (b-a) * sqrt(a + b + 1)) / ((a + b + 2) * sqrt(a * b));353 } // skewness354 355 template <class RealType, class Policy>356 BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const beta_distribution<RealType, Policy>& dist)357 {358 RealType a = dist.alpha();359 RealType b = dist.beta();360 RealType a_2 = a * a;361 RealType n = 6 * (a_2 * a - a_2 * (2 * b - 1) + b * b * (b + 1) - 2 * a * b * (b + 2));362 RealType d = a * b * (a + b + 2) * (a + b + 3);363 return n / d;364 } // kurtosis_excess365 366 template <class RealType, class Policy>367 BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const beta_distribution<RealType, Policy>& dist)368 {369 return 3 + kurtosis_excess(dist);370 } // kurtosis371 372 template <class RealType, class Policy>373 BOOST_MATH_GPU_ENABLED inline RealType pdf(const beta_distribution<RealType, Policy>& dist, const RealType& x)374 { // Probability Density/Mass Function.375 BOOST_FPU_EXCEPTION_GUARD376 377 constexpr auto function = "boost::math::pdf(beta_distribution<%1%> const&, %1%)";378 379 BOOST_MATH_STD_USING // for ADL of std functions380 381 RealType a = dist.alpha();382 RealType b = dist.beta();383 384 // Argument checks:385 RealType result = 0;386 if(false == beta_detail::check_dist_and_x(387 function,388 a, b, x,389 &result, Policy()))390 {391 return result;392 }393 using boost::math::beta;394 395 // Corner cases: check_x ensures x element of [0, 1], but PDF is 0 for x = 0 and x = 1. PDF EQN:396 // https://wikimedia.org/api/rest_v1/media/math/render/svg/125fdaa41844a8703d1a8610ac00fbf3edacc8e7397 if(x == 0)398 {399 if (a == 1)400 {401 return static_cast<RealType>(1 / beta(a, b));402 }403 else if (a < 1)404 {405 policies::raise_overflow_error<RealType>(function, nullptr, Policy());406 }407 else408 {409 return RealType(0);410 }411 }412 else if (x == 1)413 {414 if (b == 1)415 {416 return static_cast<RealType>(1 / beta(a, b));417 }418 else if (b < 1)419 {420 policies::raise_overflow_error<RealType>(function, nullptr, Policy());421 }422 else423 {424 return RealType(0);425 }426 }427 428 return static_cast<RealType>(ibeta_derivative(a, b, x, Policy()));429 } // pdf430 431 template <class RealType, class Policy>432 BOOST_MATH_GPU_ENABLED inline RealType cdf(const beta_distribution<RealType, Policy>& dist, const RealType& x)433 { // Cumulative Distribution Function beta.434 BOOST_MATH_STD_USING // for ADL of std functions435 436 constexpr auto function = "boost::math::cdf(beta_distribution<%1%> const&, %1%)";437 438 RealType a = dist.alpha();439 RealType b = dist.beta();440 441 // Argument checks:442 RealType result = 0;443 if(false == beta_detail::check_dist_and_x(444 function,445 a, b, x,446 &result, Policy()))447 {448 return result;449 }450 // Special cases:451 if (x == 0)452 {453 return 0;454 }455 else if (x == 1)456 {457 return 1;458 }459 return static_cast<RealType>(ibeta(a, b, x, Policy()));460 } // beta cdf461 462 template <class RealType, class Policy>463 BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<beta_distribution<RealType, Policy>, RealType>& c)464 { // Complemented Cumulative Distribution Function beta.465 466 BOOST_MATH_STD_USING // for ADL of std functions467 468 constexpr auto function = "boost::math::cdf(beta_distribution<%1%> const&, %1%)";469 470 RealType const& x = c.param;471 beta_distribution<RealType, Policy> const& dist = c.dist;472 RealType a = dist.alpha();473 RealType b = dist.beta();474 475 // Argument checks:476 RealType result = 0;477 if(false == beta_detail::check_dist_and_x(478 function,479 a, b, x,480 &result, Policy()))481 {482 return result;483 }484 if (x == 0)485 {486 return RealType(1);487 }488 else if (x == 1)489 {490 return RealType(0);491 }492 // Calculate cdf beta using the incomplete beta function.493 // Use of ibeta here prevents cancellation errors in calculating494 // 1 - x if x is very small, perhaps smaller than machine epsilon.495 return static_cast<RealType>(ibetac(a, b, x, Policy()));496 } // beta cdf497 498 template <class RealType, class Policy>499 BOOST_MATH_GPU_ENABLED inline RealType quantile(const beta_distribution<RealType, Policy>& dist, const RealType& p)500 { // Quantile or Percent Point beta function or501 // Inverse Cumulative probability distribution function CDF.502 // Return x (0 <= x <= 1),503 // for a given probability p (0 <= p <= 1).504 // These functions take a probability as an argument505 // and return a value such that the probability that a random variable x506 // will be less than or equal to that value507 // is whatever probability you supplied as an argument.508 509 constexpr auto function = "boost::math::quantile(beta_distribution<%1%> const&, %1%)";510 511 RealType result = 0; // of argument checks:512 RealType a = dist.alpha();513 RealType b = dist.beta();514 if(false == beta_detail::check_dist_and_prob(515 function,516 a, b, p,517 &result, Policy()))518 {519 return result;520 }521 // Special cases:522 if (p == 0)523 {524 return RealType(0);525 }526 if (p == 1)527 {528 return RealType(1);529 }530 return static_cast<RealType>(ibeta_inv(a, b, p, static_cast<RealType*>(nullptr), Policy()));531 } // quantile532 533 template <class RealType, class Policy>534 BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<beta_distribution<RealType, Policy>, RealType>& c)535 { // Complement Quantile or Percent Point beta function .536 // Return the number of expected x for a given537 // complement of the probability q.538 539 constexpr auto function = "boost::math::quantile(beta_distribution<%1%> const&, %1%)";540 541 //542 // Error checks:543 RealType q = c.param;544 const beta_distribution<RealType, Policy>& dist = c.dist;545 RealType result = 0;546 RealType a = dist.alpha();547 RealType b = dist.beta();548 if(false == beta_detail::check_dist_and_prob(549 function,550 a,551 b,552 q,553 &result, Policy()))554 {555 return result;556 }557 // Special cases:558 if(q == 1)559 {560 return RealType(0);561 }562 if(q == 0)563 {564 return RealType(1);565 }566 567 return static_cast<RealType>(ibetac_inv(a, b, q, static_cast<RealType*>(nullptr), Policy()));568 } // Quantile Complement569 570 } // namespace math571} // namespace boost572 573// This include must be at the end, *after* the accessors574// for this distribution have been defined, in order to575// keep compilers that support two-phase lookup happy.576#include <boost/math/distributions/detail/derived_accessors.hpp>577 578#if defined (BOOST_MSVC)579# pragma warning(pop)580#endif581 582#endif // BOOST_MATH_DIST_BETA_HPP583 584 585