boost/math/distributions/beta.hpp

   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>
  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>
  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>
  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>
  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>
 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>
 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           && beta_detail::check_x(function, x, result, pol);
 113       } // bool check_dist_and_x
 114
 115       template <class RealType, class Policy>
 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>
 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>
 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
 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:
 169       RealType alpha() const
 170       {
 171         return m_alpha;
 172       }
 173       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
 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
 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
 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
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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