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

source code of boost/libs/math/include/boost/math/distributions/beta.hpp