1 // boost\math\distributions\beta.hpp 2 3 // Copyright John Maddock 2006. 4 // Copyright Paul A. Bristow 2006. 5 6 // Use, modification and distribution are subject to the 7 // Boost Software License, Version 1.0. 8 // (See accompanying file LICENSE_1_0.txt 9 // or copy at http://www.boost.org/LICENSE_1_0.txt) 10 11 // http://en.wikipedia.org/wiki/Beta_distribution 12 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366h.htm 13 // http://mathworld.wolfram.com/BetaDistribution.html 14 15 // The Beta Distribution is a continuous probability distribution. 16 // The beta distribution is used to model events which are constrained to take place 17 // within an interval defined by maxima and minima, 18 // so is used extensively in PERT and other project management systems 19 // to describe the time to completion. 20 // The cdf of the beta distribution is used as a convenient way 21 // of obtaining the sum over a set of binomial outcomes. 22 // The beta distribution is also used in Bayesian statistics. 23 24 #ifndef BOOST_MATH_DIST_BETA_HPP 25 #define BOOST_MATH_DIST_BETA_HPP 26 27 #include <boost/math/distributions/fwd.hpp> 28 #include <boost/math/special_functions/beta.hpp> // for beta. 29 #include <boost/math/distributions/complement.hpp> // complements. 30 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks 31 #include <boost/math/special_functions/fpclassify.hpp> // isnan. 32 #include <boost/math/tools/roots.hpp> // for root finding. 33 34 #if defined (BOOST_MSVC) 35 # pragma warning(push) 36 # pragma warning(disable: 4702) // unreachable code 37 // in domain_error_imp in error_handling 38 #endif 39 40 #include <utility> 41 42 namespace boost 43 { 44 namespace math 45 { 46 namespace beta_detail 47 { 48 // Common error checking routines for beta distribution functions: 49 template <class RealType, class Policy> check_alpha(const char * function,const RealType & alpha,RealType * result,const Policy & pol)50 inline bool check_alpha(const char* function, const RealType& alpha, RealType* result, const Policy& pol) 51 { 52 if(!(boost::math::isfinite)(alpha) || (alpha <= 0)) 53 { 54 *result = policies::raise_domain_error<RealType>( 55 function, 56 "Alpha argument is %1%, but must be > 0 !", alpha, pol); 57 return false; 58 } 59 return true; 60 } // bool check_alpha 61 62 template <class RealType, class Policy> check_beta(const char * function,const RealType & beta,RealType * result,const Policy & pol)63 inline bool check_beta(const char* function, const RealType& beta, RealType* result, const Policy& pol) 64 { 65 if(!(boost::math::isfinite)(beta) || (beta <= 0)) 66 { 67 *result = policies::raise_domain_error<RealType>( 68 function, 69 "Beta argument is %1%, but must be > 0 !", beta, pol); 70 return false; 71 } 72 return true; 73 } // bool check_beta 74 75 template <class RealType, class Policy> check_prob(const char * function,const RealType & p,RealType * result,const Policy & pol)76 inline bool check_prob(const char* function, const RealType& p, RealType* result, const Policy& pol) 77 { 78 if((p < 0) || (p > 1) || !(boost::math::isfinite)(p)) 79 { 80 *result = policies::raise_domain_error<RealType>( 81 function, 82 "Probability argument is %1%, but must be >= 0 and <= 1 !", p, pol); 83 return false; 84 } 85 return true; 86 } // bool check_prob 87 88 template <class RealType, class Policy> check_x(const char * function,const RealType & x,RealType * result,const Policy & pol)89 inline bool check_x(const char* function, const RealType& x, RealType* result, const Policy& pol) 90 { 91 if(!(boost::math::isfinite)(x) || (x < 0) || (x > 1)) 92 { 93 *result = policies::raise_domain_error<RealType>( 94 function, 95 "x argument is %1%, but must be >= 0 and <= 1 !", x, pol); 96 return false; 97 } 98 return true; 99 } // bool check_x 100 101 template <class RealType, class Policy> check_dist(const char * function,const RealType & alpha,const RealType & beta,RealType * result,const Policy & pol)102 inline bool check_dist(const char* function, const RealType& alpha, const RealType& beta, RealType* result, const Policy& pol) 103 { // Check both alpha and beta. 104 return check_alpha(function, alpha, result, pol) 105 && check_beta(function, beta, result, pol); 106 } // bool check_dist 107 108 template <class RealType, class Policy> check_dist_and_x(const char * function,const RealType & alpha,const RealType & beta,RealType x,RealType * result,const Policy & pol)109 inline bool check_dist_and_x(const char* function, const RealType& alpha, const RealType& beta, RealType x, RealType* result, const Policy& pol) 110 { 111 return check_dist(function, alpha, beta, result, pol) 112 && check_x(function, x, result, pol); 113 } // bool check_dist_and_x 114 115 template <class RealType, class Policy> check_dist_and_prob(const char * function,const RealType & alpha,const RealType & beta,RealType p,RealType * result,const Policy & pol)116 inline bool check_dist_and_prob(const char* function, const RealType& alpha, const RealType& beta, RealType p, RealType* result, const Policy& pol) 117 { 118 return check_dist(function, alpha, beta, result, pol) 119 && check_prob(function, p, result, pol); 120 } // bool check_dist_and_prob 121 122 template <class RealType, class Policy> check_mean(const char * function,const RealType & mean,RealType * result,const Policy & pol)123 inline bool check_mean(const char* function, const RealType& mean, RealType* result, const Policy& pol) 124 { 125 if(!(boost::math::isfinite)(mean) || (mean <= 0)) 126 { 127 *result = policies::raise_domain_error<RealType>( 128 function, 129 "mean argument is %1%, but must be > 0 !", mean, pol); 130 return false; 131 } 132 return true; 133 } // bool check_mean 134 template <class RealType, class Policy> check_variance(const char * function,const RealType & variance,RealType * result,const Policy & pol)135 inline bool check_variance(const char* function, const RealType& variance, RealType* result, const Policy& pol) 136 { 137 if(!(boost::math::isfinite)(variance) || (variance <= 0)) 138 { 139 *result = policies::raise_domain_error<RealType>( 140 function, 141 "variance argument is %1%, but must be > 0 !", variance, pol); 142 return false; 143 } 144 return true; 145 } // bool check_variance 146 } // namespace beta_detail 147 148 // typedef beta_distribution<double> beta; 149 // is deliberately NOT included to avoid a name clash with the beta function. 150 // Use beta_distribution<> mybeta(...) to construct type double. 151 152 template <class RealType = double, class Policy = policies::policy<> > 153 class beta_distribution 154 { 155 public: 156 typedef RealType value_type; 157 typedef Policy policy_type; 158 beta_distribution(RealType alpha=1,RealType beta=1)159 beta_distribution(RealType alpha = 1, RealType beta = 1) : m_alpha(alpha), m_beta(beta) 160 { 161 RealType result; 162 beta_detail::check_dist( 163 "boost::math::beta_distribution<%1%>::beta_distribution", 164 m_alpha, 165 m_beta, 166 &result, Policy()); 167 } // beta_distribution constructor. 168 // Accessor functions: alpha() const169 RealType alpha() const 170 { 171 return m_alpha; 172 } beta() const173 RealType beta() const 174 { // . 175 return m_beta; 176 } 177 178 // Estimation of the alpha & beta parameters. 179 // http://en.wikipedia.org/wiki/Beta_distribution 180 // gives formulae in section on parameter estimation. 181 // Also NIST EDA page 3 & 4 give the same. 182 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366h.htm 183 // http://www.epi.ucdavis.edu/diagnostictests/betabuster.html 184 find_alpha(RealType mean,RealType variance)185 static RealType find_alpha( 186 RealType mean, // Expected value of mean. 187 RealType variance) // Expected value of variance. 188 { 189 static const char* function = "boost::math::beta_distribution<%1%>::find_alpha"; 190 RealType result = 0; // of error checks. 191 if(false == 192 beta_detail::check_mean( 193 function, mean, &result, Policy()) 194 && 195 beta_detail::check_variance( 196 function, variance, &result, Policy()) 197 ) 198 { 199 return result; 200 } 201 return mean * (( (mean * (1 - mean)) / variance)- 1); 202 } // RealType find_alpha 203 find_beta(RealType mean,RealType variance)204 static RealType find_beta( 205 RealType mean, // Expected value of mean. 206 RealType variance) // Expected value of variance. 207 { 208 static const char* function = "boost::math::beta_distribution<%1%>::find_beta"; 209 RealType result = 0; // of error checks. 210 if(false == 211 beta_detail::check_mean( 212 function, mean, &result, Policy()) 213 && 214 beta_detail::check_variance( 215 function, variance, &result, Policy()) 216 ) 217 { 218 return result; 219 } 220 return (1 - mean) * (((mean * (1 - mean)) /variance)-1); 221 } // RealType find_beta 222 223 // Estimate alpha & beta from either alpha or beta, and x and probability. 224 // Uses for these parameter estimators are unclear. 225 find_alpha(RealType beta,RealType x,RealType probability)226 static RealType find_alpha( 227 RealType beta, // from beta. 228 RealType x, // x. 229 RealType probability) // cdf 230 { 231 static const char* function = "boost::math::beta_distribution<%1%>::find_alpha"; 232 RealType result = 0; // of error checks. 233 if(false == 234 beta_detail::check_prob( 235 function, probability, &result, Policy()) 236 && 237 beta_detail::check_beta( 238 function, beta, &result, Policy()) 239 && 240 beta_detail::check_x( 241 function, x, &result, Policy()) 242 ) 243 { 244 return result; 245 } 246 return ibeta_inva(beta, x, probability, Policy()); 247 } // RealType find_alpha(beta, a, probability) 248 find_beta(RealType alpha,RealType x,RealType probability)249 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 static const char* function = "boost::math::beta_distribution<%1%>::find_beta"; 256 RealType result = 0; // of error checks. 257 if(false == 258 beta_detail::check_prob( 259 function, probability, &result, Policy()) 260 && 261 beta_detail::check_alpha( 262 function, alpha, &result, Policy()) 263 && 264 beta_detail::check_x( 265 function, x, &result, Policy()) 266 ) 267 { 268 return result; 269 } 270 return ibeta_invb(alpha, x, probability, Policy()); 271 } // RealType find_beta(alpha, x, probability) 272 273 private: 274 RealType m_alpha; // Two parameters of the beta distribution. 275 RealType m_beta; 276 }; // template <class RealType, class Policy> class beta_distribution 277 278 template <class RealType, class Policy> range(const beta_distribution<RealType,Policy> &)279 inline const std::pair<RealType, RealType> range(const beta_distribution<RealType, Policy>& /* dist */) 280 { // Range of permissible values for random variable x. 281 using boost::math::tools::max_value; 282 return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1)); 283 } 284 285 template <class RealType, class Policy> support(const beta_distribution<RealType,Policy> &)286 inline const std::pair<RealType, RealType> support(const beta_distribution<RealType, Policy>& /* dist */) 287 { // Range of supported values for random variable x. 288 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. 289 return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1)); 290 } 291 292 template <class RealType, class Policy> mean(const beta_distribution<RealType,Policy> & dist)293 inline RealType mean(const beta_distribution<RealType, Policy>& dist) 294 { // Mean of beta distribution = np. 295 return dist.alpha() / (dist.alpha() + dist.beta()); 296 } // mean 297 298 template <class RealType, class Policy> variance(const beta_distribution<RealType,Policy> & dist)299 inline RealType variance(const beta_distribution<RealType, Policy>& dist) 300 { // Variance of beta distribution = np(1-p). 301 RealType a = dist.alpha(); 302 RealType b = dist.beta(); 303 return (a * b) / ((a + b ) * (a + b) * (a + b + 1)); 304 } // variance 305 306 template <class RealType, class Policy> mode(const beta_distribution<RealType,Policy> & dist)307 inline RealType mode(const beta_distribution<RealType, Policy>& dist) 308 { 309 static const char* function = "boost::math::mode(beta_distribution<%1%> const&)"; 310 311 RealType result; 312 if ((dist.alpha() <= 1)) 313 { 314 result = policies::raise_domain_error<RealType>( 315 function, 316 "mode undefined for alpha = %1%, must be > 1!", dist.alpha(), Policy()); 317 return result; 318 } 319 320 if ((dist.beta() <= 1)) 321 { 322 result = policies::raise_domain_error<RealType>( 323 function, 324 "mode undefined for beta = %1%, must be > 1!", dist.beta(), Policy()); 325 return result; 326 } 327 RealType a = dist.alpha(); 328 RealType b = dist.beta(); 329 return (a-1) / (a + b - 2); 330 } // mode 331 332 //template <class RealType, class Policy> 333 //inline RealType median(const beta_distribution<RealType, Policy>& dist) 334 //{ // Median of beta distribution is not defined. 335 // return tools::domain_error<RealType>(function, "Median is not implemented, result is %1%!", std::numeric_limits<RealType>::quiet_NaN()); 336 //} // median 337 338 //But WILL be provided by the derived accessor as quantile(0.5). 339 340 template <class RealType, class Policy> skewness(const beta_distribution<RealType,Policy> & dist)341 inline RealType skewness(const beta_distribution<RealType, Policy>& dist) 342 { 343 BOOST_MATH_STD_USING // ADL of std functions. 344 RealType a = dist.alpha(); 345 RealType b = dist.beta(); 346 return (2 * (b-a) * sqrt(a + b + 1)) / ((a + b + 2) * sqrt(a * b)); 347 } // skewness 348 349 template <class RealType, class Policy> kurtosis_excess(const beta_distribution<RealType,Policy> & dist)350 inline RealType kurtosis_excess(const beta_distribution<RealType, Policy>& dist) 351 { 352 RealType a = dist.alpha(); 353 RealType b = dist.beta(); 354 RealType a_2 = a * a; 355 RealType n = 6 * (a_2 * a - a_2 * (2 * b - 1) + b * b * (b + 1) - 2 * a * b * (b + 2)); 356 RealType d = a * b * (a + b + 2) * (a + b + 3); 357 return n / d; 358 } // kurtosis_excess 359 360 template <class RealType, class Policy> kurtosis(const beta_distribution<RealType,Policy> & dist)361 inline RealType kurtosis(const beta_distribution<RealType, Policy>& dist) 362 { 363 return 3 + kurtosis_excess(dist); 364 } // kurtosis 365 366 template <class RealType, class Policy> pdf(const beta_distribution<RealType,Policy> & dist,const RealType & x)367 inline RealType pdf(const beta_distribution<RealType, Policy>& dist, const RealType& x) 368 { // Probability Density/Mass Function. 369 BOOST_FPU_EXCEPTION_GUARD 370 371 static const char* function = "boost::math::pdf(beta_distribution<%1%> const&, %1%)"; 372 373 BOOST_MATH_STD_USING // for ADL of std functions 374 375 RealType a = dist.alpha(); 376 RealType b = dist.beta(); 377 378 // Argument checks: 379 RealType result = 0; 380 if(false == beta_detail::check_dist_and_x( 381 function, 382 a, b, x, 383 &result, Policy())) 384 { 385 return result; 386 } 387 using boost::math::beta; 388 return ibeta_derivative(a, b, x, Policy()); 389 } // pdf 390 391 template <class RealType, class Policy> cdf(const beta_distribution<RealType,Policy> & dist,const RealType & x)392 inline RealType cdf(const beta_distribution<RealType, Policy>& dist, const RealType& x) 393 { // Cumulative Distribution Function beta. 394 BOOST_MATH_STD_USING // for ADL of std functions 395 396 static const char* function = "boost::math::cdf(beta_distribution<%1%> const&, %1%)"; 397 398 RealType a = dist.alpha(); 399 RealType b = dist.beta(); 400 401 // Argument checks: 402 RealType result = 0; 403 if(false == beta_detail::check_dist_and_x( 404 function, 405 a, b, x, 406 &result, Policy())) 407 { 408 return result; 409 } 410 // Special cases: 411 if (x == 0) 412 { 413 return 0; 414 } 415 else if (x == 1) 416 { 417 return 1; 418 } 419 return ibeta(a, b, x, Policy()); 420 } // beta cdf 421 422 template <class RealType, class Policy> cdf(const complemented2_type<beta_distribution<RealType,Policy>,RealType> & c)423 inline RealType cdf(const complemented2_type<beta_distribution<RealType, Policy>, RealType>& c) 424 { // Complemented Cumulative Distribution Function beta. 425 426 BOOST_MATH_STD_USING // for ADL of std functions 427 428 static const char* function = "boost::math::cdf(beta_distribution<%1%> const&, %1%)"; 429 430 RealType const& x = c.param; 431 beta_distribution<RealType, Policy> const& dist = c.dist; 432 RealType a = dist.alpha(); 433 RealType b = dist.beta(); 434 435 // Argument checks: 436 RealType result = 0; 437 if(false == beta_detail::check_dist_and_x( 438 function, 439 a, b, x, 440 &result, Policy())) 441 { 442 return result; 443 } 444 if (x == 0) 445 { 446 return 1; 447 } 448 else if (x == 1) 449 { 450 return 0; 451 } 452 // Calculate cdf beta using the incomplete beta function. 453 // Use of ibeta here prevents cancellation errors in calculating 454 // 1 - x if x is very small, perhaps smaller than machine epsilon. 455 return ibetac(a, b, x, Policy()); 456 } // beta cdf 457 458 template <class RealType, class Policy> quantile(const beta_distribution<RealType,Policy> & dist,const RealType & p)459 inline RealType quantile(const beta_distribution<RealType, Policy>& dist, const RealType& p) 460 { // Quantile or Percent Point beta function or 461 // Inverse Cumulative probability distribution function CDF. 462 // Return x (0 <= x <= 1), 463 // for a given probability p (0 <= p <= 1). 464 // These functions take a probability as an argument 465 // and return a value such that the probability that a random variable x 466 // will be less than or equal to that value 467 // is whatever probability you supplied as an argument. 468 469 static const char* function = "boost::math::quantile(beta_distribution<%1%> const&, %1%)"; 470 471 RealType result = 0; // of argument checks: 472 RealType a = dist.alpha(); 473 RealType b = dist.beta(); 474 if(false == beta_detail::check_dist_and_prob( 475 function, 476 a, b, p, 477 &result, Policy())) 478 { 479 return result; 480 } 481 // Special cases: 482 if (p == 0) 483 { 484 return 0; 485 } 486 if (p == 1) 487 { 488 return 1; 489 } 490 return ibeta_inv(a, b, p, static_cast<RealType*>(0), Policy()); 491 } // quantile 492 493 template <class RealType, class Policy> quantile(const complemented2_type<beta_distribution<RealType,Policy>,RealType> & c)494 inline RealType quantile(const complemented2_type<beta_distribution<RealType, Policy>, RealType>& c) 495 { // Complement Quantile or Percent Point beta function . 496 // Return the number of expected x for a given 497 // complement of the probability q. 498 499 static const char* function = "boost::math::quantile(beta_distribution<%1%> const&, %1%)"; 500 501 // 502 // Error checks: 503 RealType q = c.param; 504 const beta_distribution<RealType, Policy>& dist = c.dist; 505 RealType result = 0; 506 RealType a = dist.alpha(); 507 RealType b = dist.beta(); 508 if(false == beta_detail::check_dist_and_prob( 509 function, 510 a, 511 b, 512 q, 513 &result, Policy())) 514 { 515 return result; 516 } 517 // Special cases: 518 if(q == 1) 519 { 520 return 0; 521 } 522 if(q == 0) 523 { 524 return 1; 525 } 526 527 return ibetac_inv(a, b, q, static_cast<RealType*>(0), Policy()); 528 } // Quantile Complement 529 530 } // namespace math 531 } // namespace boost 532 533 // This include must be at the end, *after* the accessors 534 // for this distribution have been defined, in order to 535 // keep compilers that support two-phase lookup happy. 536 #include <boost/math/distributions/detail/derived_accessors.hpp> 537 538 #if defined (BOOST_MSVC) 539 # pragma warning(pop) 540 #endif 541 542 #endif // BOOST_MATH_DIST_BETA_HPP 543 544 545