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