math/distributions/beta.hpp

0001 // boost\math\distributions\beta.hpp
0002
0003 // Copyright John Maddock 2006.
0004 // Copyright Paul A. Bristow 2006.
0005 // Copyright Matt Borland 2023.
0006
0007 // Use, modification and distribution are subject to the
0008 // Boost Software License, Version 1.0.
0009 // (See accompanying file LICENSE_1_0.txt
0010 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0011
0012 // http://en.wikipedia.org/wiki/Beta_distribution
0013 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366h.htm
0014 // http://mathworld.wolfram.com/BetaDistribution.html
0015
0016 // The Beta Distribution is a continuous probability distribution.
0017 // The beta distribution is used to model events which are constrained to take place
0018 // within an interval defined by maxima and minima,
0019 // so is used extensively in PERT and other project management systems
0020 // to describe the time to completion.
0021 // The cdf of the beta distribution is used as a convenient way
0022 // of obtaining the sum over a set of binomial outcomes.
0023 // The beta distribution is also used in Bayesian statistics.
0024
0025 #ifndef BOOST_MATH_DIST_BETA_HPP
0026 #define BOOST_MATH_DIST_BETA_HPP
0027
0028 #include <boost/math/distributions/fwd.hpp>
0029 #include <boost/math/special_functions/beta.hpp> // for beta.
0030 #include <boost/math/distributions/complement.hpp> // complements.
0031 #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
0032 #include <boost/math/special_functions/fpclassify.hpp> // isnan.
0033 #include <boost/math/tools/roots.hpp> // for root finding.
0034
0035 #if defined (BOOST_MSVC)
0036 #  pragma warning(push)
0037 #  pragma warning(disable: 4702) // unreachable code
0038 // in domain_error_imp in error_handling
0039 #endif
0040
0041 #include <utility>
0042
0043 namespace boost
0044 {
0045   namespace math
0046   {
0047     namespace beta_detail
0048     {
0049       // Common error checking routines for beta distribution functions:
0050       template <class RealType, class Policy>
0051       inline bool check_alpha(const char* function, const RealType& alpha, RealType* result, const Policy& pol)
0052       {
0053         if(!(boost::math::isfinite)(alpha) || (alpha <= 0))
0054         {
0055           *result = policies::raise_domain_error<RealType>(
0056             function,
0057             "Alpha argument is %1%, but must be > 0 !", alpha, pol);
0058           return false;
0059         }
0060         return true;
0061       } // bool check_alpha
0062
0063       template <class RealType, class Policy>
0064       inline bool check_beta(const char* function, const RealType& beta, RealType* result, const Policy& pol)
0065       {
0066         if(!(boost::math::isfinite)(beta) || (beta <= 0))
0067         {
0068           *result = policies::raise_domain_error<RealType>(
0069             function,
0070             "Beta argument is %1%, but must be > 0 !", beta, pol);
0071           return false;
0072         }
0073         return true;
0074       } // bool check_beta
0075
0076       template <class RealType, class Policy>
0077       inline bool check_prob(const char* function, const RealType& p, RealType* result, const Policy& pol)
0078       {
0079         if((p < 0) || (p > 1) || !(boost::math::isfinite)(p))
0080         {
0081           *result = policies::raise_domain_error<RealType>(
0082             function,
0083             "Probability argument is %1%, but must be >= 0 and <= 1 !", p, pol);
0084           return false;
0085         }
0086         return true;
0087       } // bool check_prob
0088
0089       template <class RealType, class Policy>
0090       inline bool check_x(const char* function, const RealType& x, RealType* result, const Policy& pol)
0091       {
0092         if(!(boost::math::isfinite)(x) || (x < 0) || (x > 1))
0093         {
0094           *result = policies::raise_domain_error<RealType>(
0095             function,
0096             "x argument is %1%, but must be >= 0 and <= 1 !", x, pol);
0097           return false;
0098         }
0099         return true;
0100       } // bool check_x
0101
0102       template <class RealType, class Policy>
0103       inline bool check_dist(const char* function, const RealType& alpha, const RealType& beta, RealType* result, const Policy& pol)
0104       { // Check both alpha and beta.
0105         return check_alpha(function, alpha, result, pol)
0106           && check_beta(function, beta, result, pol);
0107       } // bool check_dist
0108
0109       template <class RealType, class Policy>
0110       inline bool check_dist_and_x(const char* function, const RealType& alpha, const RealType& beta, RealType x, RealType* result, const Policy& pol)
0111       {
0112         return check_dist(function, alpha, beta, result, pol)
0113           && beta_detail::check_x(function, x, result, pol);
0114       } // bool check_dist_and_x
0115
0116       template <class RealType, class Policy>
0117       inline bool check_dist_and_prob(const char* function, const RealType& alpha, const RealType& beta, RealType p, RealType* result, const Policy& pol)
0118       {
0119         return check_dist(function, alpha, beta, result, pol)
0120           && check_prob(function, p, result, pol);
0121       } // bool check_dist_and_prob
0122
0123       template <class RealType, class Policy>
0124       inline bool check_mean(const char* function, const RealType& mean, RealType* result, const Policy& pol)
0125       {
0126         if(!(boost::math::isfinite)(mean) || (mean <= 0))
0127         {
0128           *result = policies::raise_domain_error<RealType>(
0129             function,
0130             "mean argument is %1%, but must be > 0 !", mean, pol);
0131           return false;
0132         }
0133         return true;
0134       } // bool check_mean
0135       template <class RealType, class Policy>
0136       inline bool check_variance(const char* function, const RealType& variance, RealType* result, const Policy& pol)
0137       {
0138         if(!(boost::math::isfinite)(variance) || (variance <= 0))
0139         {
0140           *result = policies::raise_domain_error<RealType>(
0141             function,
0142             "variance argument is %1%, but must be > 0 !", variance, pol);
0143           return false;
0144         }
0145         return true;
0146       } // bool check_variance
0147     } // namespace beta_detail
0148
0149     // typedef beta_distribution<double> beta;
0150     // is deliberately NOT included to avoid a name clash with the beta function.
0151     // Use beta_distribution<> mybeta(...) to construct type double.
0152
0153     template <class RealType = double, class Policy = policies::policy<> >
0154     class beta_distribution
0155     {
0156     public:
0157       typedef RealType value_type;
0158       typedef Policy policy_type;
0159
0160       beta_distribution(RealType l_alpha = 1, RealType l_beta = 1) : m_alpha(l_alpha), m_beta(l_beta)
0161       {
0162         RealType result;
0163         beta_detail::check_dist(
0164            "boost::math::beta_distribution<%1%>::beta_distribution",
0165           m_alpha,
0166           m_beta,
0167           &result, Policy());
0168       } // beta_distribution constructor.
0169       // Accessor functions:
0170       RealType alpha() const
0171       {
0172         return m_alpha;
0173       }
0174       RealType beta() const
0175       { // .
0176         return m_beta;
0177       }
0178
0179       // Estimation of the alpha & beta parameters.
0180       // http://en.wikipedia.org/wiki/Beta_distribution
0181       // gives formulae in section on parameter estimation.
0182       // Also NIST EDA page 3 & 4 give the same.
0183       // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366h.htm
0184       // http://www.epi.ucdavis.edu/diagnostictests/betabuster.html
0185
0186       static RealType find_alpha(
0187         RealType mean, // Expected value of mean.
0188         RealType variance) // Expected value of variance.
0189       {
0190         static const char* function = "boost::math::beta_distribution<%1%>::find_alpha";
0191         RealType result = 0; // of error checks.
0192         if(false ==
0193             (
0194               beta_detail::check_mean(function, mean, &result, Policy())
0195               && beta_detail::check_variance(function, variance, &result, Policy())
0196             )
0197           )
0198         {
0199           return result;
0200         }
0201         return mean * (( (mean * (1 - mean)) / variance)- 1);
0202       } // RealType find_alpha
0203
0204       static RealType find_beta(
0205         RealType mean, // Expected value of mean.
0206         RealType variance) // Expected value of variance.
0207       {
0208         static const char* function = "boost::math::beta_distribution<%1%>::find_beta";
0209         RealType result = 0; // of error checks.
0210         if(false ==
0211             (
0212               beta_detail::check_mean(function, mean, &result, Policy())
0213               &&
0214               beta_detail::check_variance(function, variance, &result, Policy())
0215             )
0216           )
0217         {
0218           return result;
0219         }
0220         return (1 - mean) * (((mean * (1 - mean)) /variance)-1);
0221       } //  RealType find_beta
0222
0223       // Estimate alpha & beta from either alpha or beta, and x and probability.
0224       // Uses for these parameter estimators are unclear.
0225
0226       static RealType find_alpha(
0227         RealType beta, // from beta.
0228         RealType x, //  x.
0229         RealType probability) // cdf
0230       {
0231         static const char* function = "boost::math::beta_distribution<%1%>::find_alpha";
0232         RealType result = 0; // of error checks.
0233         if(false ==
0234             (
0235              beta_detail::check_prob(function, probability, &result, Policy())
0236              &&
0237              beta_detail::check_beta(function, beta, &result, Policy())
0238              &&
0239              beta_detail::check_x(function, x, &result, Policy())
0240             )
0241           )
0242         {
0243           return result;
0244         }
0245         return static_cast<RealType>(ibeta_inva(beta, x, probability, Policy()));
0246       } // RealType find_alpha(beta, a, probability)
0247
0248       static RealType find_beta(
0249         // ibeta_invb(T b, T x, T p); (alpha, x, cdf,)
0250         RealType alpha, // alpha.
0251         RealType x, // probability x.
0252         RealType probability) // probability cdf.
0253       {
0254         static const char* function = "boost::math::beta_distribution<%1%>::find_beta";
0255         RealType result = 0; // of error checks.
0256         if(false ==
0257             (
0258               beta_detail::check_prob(function, probability, &result, Policy())
0259               &&
0260               beta_detail::check_alpha(function, alpha, &result, Policy())
0261               &&
0262               beta_detail::check_x(function, x, &result, Policy())
0263             )
0264           )
0265         {
0266           return result;
0267         }
0268         return static_cast<RealType>(ibeta_invb(alpha, x, probability, Policy()));
0269       } //  RealType find_beta(alpha, x, probability)
0270
0271     private:
0272       RealType m_alpha; // Two parameters of the beta distribution.
0273       RealType m_beta;
0274     }; // template <class RealType, class Policy> class beta_distribution
0275
0276     #ifdef __cpp_deduction_guides
0277     template <class RealType>
0278     beta_distribution(RealType)->beta_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0279     template <class RealType>
0280     beta_distribution(RealType, RealType)->beta_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0281     #endif
0282
0283     template <class RealType, class Policy>
0284     inline const std::pair<RealType, RealType> range(const beta_distribution<RealType, Policy>& /* dist */)
0285     { // Range of permissible values for random variable x.
0286       using boost::math::tools::max_value;
0287       return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1));
0288     }
0289
0290     template <class RealType, class Policy>
0291     inline const std::pair<RealType, RealType> support(const beta_distribution<RealType, Policy>&  /* dist */)
0292     { // Range of supported values for random variable x.
0293       // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0294       return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1));
0295     }
0296
0297     template <class RealType, class Policy>
0298     inline RealType mean(const beta_distribution<RealType, Policy>& dist)
0299     { // Mean of beta distribution = np.
0300       return  dist.alpha() / (dist.alpha() + dist.beta());
0301     } // mean
0302
0303     template <class RealType, class Policy>
0304     inline RealType variance(const beta_distribution<RealType, Policy>& dist)
0305     { // Variance of beta distribution = np(1-p).
0306       RealType a = dist.alpha();
0307       RealType b = dist.beta();
0308       return  (a * b) / ((a + b ) * (a + b) * (a + b + 1));
0309     } // variance
0310
0311     template <class RealType, class Policy>
0312     inline RealType mode(const beta_distribution<RealType, Policy>& dist)
0313     {
0314       static const char* function = "boost::math::mode(beta_distribution<%1%> const&)";
0315
0316       RealType result;
0317       if ((dist.alpha() <= 1))
0318       {
0319         result = policies::raise_domain_error<RealType>(
0320           function,
0321           "mode undefined for alpha = %1%, must be > 1!", dist.alpha(), Policy());
0322         return result;
0323       }
0324
0325       if ((dist.beta() <= 1))
0326       {
0327         result = policies::raise_domain_error<RealType>(
0328           function,
0329           "mode undefined for beta = %1%, must be > 1!", dist.beta(), Policy());
0330         return result;
0331       }
0332       RealType a = dist.alpha();
0333       RealType b = dist.beta();
0334       return (a-1) / (a + b - 2);
0335     } // mode
0336
0337     //template <class RealType, class Policy>
0338     //inline RealType median(const beta_distribution<RealType, Policy>& dist)
0339     //{ // Median of beta distribution is not defined.
0340     //  return tools::domain_error<RealType>(function, "Median is not implemented, result is %1%!", std::numeric_limits<RealType>::quiet_NaN());
0341     //} // median
0342
0343     //But WILL be provided by the derived accessor as quantile(0.5).
0344
0345     template <class RealType, class Policy>
0346     inline RealType skewness(const beta_distribution<RealType, Policy>& dist)
0347     {
0348       BOOST_MATH_STD_USING // ADL of std functions.
0349       RealType a = dist.alpha();
0350       RealType b = dist.beta();
0351       return (2 * (b-a) * sqrt(a + b + 1)) / ((a + b + 2) * sqrt(a * b));
0352     } // skewness
0353
0354     template <class RealType, class Policy>
0355     inline RealType kurtosis_excess(const beta_distribution<RealType, Policy>& dist)
0356     {
0357       RealType a = dist.alpha();
0358       RealType b = dist.beta();
0359       RealType a_2 = a * a;
0360       RealType n = 6 * (a_2 * a - a_2 * (2 * b - 1) + b * b * (b + 1) - 2 * a * b * (b + 2));
0361       RealType d = a * b * (a + b + 2) * (a + b + 3);
0362       return  n / d;
0363     } // kurtosis_excess
0364
0365     template <class RealType, class Policy>
0366     inline RealType kurtosis(const beta_distribution<RealType, Policy>& dist)
0367     {
0368       return 3 + kurtosis_excess(dist);
0369     } // kurtosis
0370
0371     template <class RealType, class Policy>
0372     inline RealType pdf(const beta_distribution<RealType, Policy>& dist, const RealType& x)
0373     { // Probability Density/Mass Function.
0374       BOOST_FPU_EXCEPTION_GUARD
0375
0376       static const char* function = "boost::math::pdf(beta_distribution<%1%> const&, %1%)";
0377
0378       BOOST_MATH_STD_USING // for ADL of std functions
0379
0380       RealType a = dist.alpha();
0381       RealType b = dist.beta();
0382
0383       // Argument checks:
0384       RealType result = 0;
0385       if(false == beta_detail::check_dist_and_x(
0386         function,
0387         a, b, x,
0388         &result, Policy()))
0389       {
0390         return result;
0391       }
0392       using boost::math::beta;
0393
0394       // Corner cases: check_x ensures x element of [0, 1], but PDF is 0 for x = 0 and x = 1. PDF EQN:
0395       // https://wikimedia.org/api/rest_v1/media/math/render/svg/125fdaa41844a8703d1a8610ac00fbf3edacc8e7
0396       if(x == 0)
0397       {
0398         if (a == 1)
0399         {
0400           return static_cast<RealType>(1 / beta(a, b));
0401         }
0402         else if (a < 1)
0403         {
0404           policies::raise_overflow_error<RealType>(function, nullptr, Policy());
0405         }
0406         else
0407         {
0408           return RealType(0);
0409         }
0410       }
0411       else if (x == 1)
0412       {
0413         if (b == 1)
0414         {
0415           return static_cast<RealType>(1 / beta(a, b));
0416         }
0417         else if (b < 1)
0418         {
0419           policies::raise_overflow_error<RealType>(function, nullptr, Policy());
0420         }
0421         else
0422         {
0423           return RealType(0);
0424         }
0425       }
0426
0427       return static_cast<RealType>(ibeta_derivative(a, b, x, Policy()));
0428     } // pdf
0429
0430     template <class RealType, class Policy>
0431     inline RealType cdf(const beta_distribution<RealType, Policy>& dist, const RealType& x)
0432     { // Cumulative Distribution Function beta.
0433       BOOST_MATH_STD_USING // for ADL of std functions
0434
0435       static const char* function = "boost::math::cdf(beta_distribution<%1%> const&, %1%)";
0436
0437       RealType a = dist.alpha();
0438       RealType b = dist.beta();
0439
0440       // Argument checks:
0441       RealType result = 0;
0442       if(false == beta_detail::check_dist_and_x(
0443         function,
0444         a, b, x,
0445         &result, Policy()))
0446       {
0447         return result;
0448       }
0449       // Special cases:
0450       if (x == 0)
0451       {
0452         return 0;
0453       }
0454       else if (x == 1)
0455       {
0456         return 1;
0457       }
0458       return static_cast<RealType>(ibeta(a, b, x, Policy()));
0459     } // beta cdf
0460
0461     template <class RealType, class Policy>
0462     inline RealType cdf(const complemented2_type<beta_distribution<RealType, Policy>, RealType>& c)
0463     { // Complemented Cumulative Distribution Function beta.
0464
0465       BOOST_MATH_STD_USING // for ADL of std functions
0466
0467       static const char* function = "boost::math::cdf(beta_distribution<%1%> const&, %1%)";
0468
0469       RealType const& x = c.param;
0470       beta_distribution<RealType, Policy> const& dist = c.dist;
0471       RealType a = dist.alpha();
0472       RealType b = dist.beta();
0473
0474       // Argument checks:
0475       RealType result = 0;
0476       if(false == beta_detail::check_dist_and_x(
0477         function,
0478         a, b, x,
0479         &result, Policy()))
0480       {
0481         return result;
0482       }
0483       if (x == 0)
0484       {
0485         return RealType(1);
0486       }
0487       else if (x == 1)
0488       {
0489         return RealType(0);
0490       }
0491       // Calculate cdf beta using the incomplete beta function.
0492       // Use of ibeta here prevents cancellation errors in calculating
0493       // 1 - x if x is very small, perhaps smaller than machine epsilon.
0494       return static_cast<RealType>(ibetac(a, b, x, Policy()));
0495     } // beta cdf
0496
0497     template <class RealType, class Policy>
0498     inline RealType quantile(const beta_distribution<RealType, Policy>& dist, const RealType& p)
0499     { // Quantile or Percent Point beta function or
0500       // Inverse Cumulative probability distribution function CDF.
0501       // Return x (0 <= x <= 1),
0502       // for a given probability p (0 <= p <= 1).
0503       // These functions take a probability as an argument
0504       // and return a value such that the probability that a random variable x
0505       // will be less than or equal to that value
0506       // is whatever probability you supplied as an argument.
0507
0508       static const char* function = "boost::math::quantile(beta_distribution<%1%> const&, %1%)";
0509
0510       RealType result = 0; // of argument checks:
0511       RealType a = dist.alpha();
0512       RealType b = dist.beta();
0513       if(false == beta_detail::check_dist_and_prob(
0514         function,
0515         a, b, p,
0516         &result, Policy()))
0517       {
0518         return result;
0519       }
0520       // Special cases:
0521       if (p == 0)
0522       {
0523         return RealType(0);
0524       }
0525       if (p == 1)
0526       {
0527         return RealType(1);
0528       }
0529       return static_cast<RealType>(ibeta_inv(a, b, p, static_cast<RealType*>(nullptr), Policy()));
0530     } // quantile
0531
0532     template <class RealType, class Policy>
0533     inline RealType quantile(const complemented2_type<beta_distribution<RealType, Policy>, RealType>& c)
0534     { // Complement Quantile or Percent Point beta function .
0535       // Return the number of expected x for a given
0536       // complement of the probability q.
0537
0538       static const char* function = "boost::math::quantile(beta_distribution<%1%> const&, %1%)";
0539
0540       //
0541       // Error checks:
0542       RealType q = c.param;
0543       const beta_distribution<RealType, Policy>& dist = c.dist;
0544       RealType result = 0;
0545       RealType a = dist.alpha();
0546       RealType b = dist.beta();
0547       if(false == beta_detail::check_dist_and_prob(
0548         function,
0549         a,
0550         b,
0551         q,
0552         &result, Policy()))
0553       {
0554         return result;
0555       }
0556       // Special cases:
0557       if(q == 1)
0558       {
0559         return RealType(0);
0560       }
0561       if(q == 0)
0562       {
0563         return RealType(1);
0564       }
0565
0566       return static_cast<RealType>(ibetac_inv(a, b, q, static_cast<RealType*>(nullptr), Policy()));
0567     } // Quantile Complement
0568
0569   } // namespace math
0570 } // namespace boost
0571
0572 // This include must be at the end, *after* the accessors
0573 // for this distribution have been defined, in order to
0574 // keep compilers that support two-phase lookup happy.
0575 #include <boost/math/distributions/detail/derived_accessors.hpp>
0576
0577 #if defined (BOOST_MSVC)
0578 # pragma warning(pop)
0579 #endif
0580
0581 #endif // BOOST_MATH_DIST_BETA_HPP
0582
0583