]>
Commit | Line | Data |
---|---|---|
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 l_alpha = 1, RealType l_beta = 1) : m_alpha(l_alpha), m_beta(l_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 | ( | |
193 | beta_detail::check_mean(function, mean, &result, Policy()) | |
194 | && beta_detail::check_variance(function, variance, &result, Policy()) | |
195 | ) | |
196 | ) | |
197 | { | |
198 | return result; | |
199 | } | |
200 | return mean * (( (mean * (1 - mean)) / variance)- 1); | |
201 | } // RealType find_alpha | |
202 | ||
203 | static RealType find_beta( | |
204 | RealType mean, // Expected value of mean. | |
205 | RealType variance) // Expected value of variance. | |
206 | { | |
207 | static const char* function = "boost::math::beta_distribution<%1%>::find_beta"; | |
208 | RealType result = 0; // of error checks. | |
209 | if(false == | |
210 | ( | |
211 | beta_detail::check_mean(function, mean, &result, Policy()) | |
212 | && | |
213 | beta_detail::check_variance(function, variance, &result, Policy()) | |
214 | ) | |
215 | ) | |
216 | { | |
217 | return result; | |
218 | } | |
219 | return (1 - mean) * (((mean * (1 - mean)) /variance)-1); | |
220 | } // RealType find_beta | |
221 | ||
222 | // Estimate alpha & beta from either alpha or beta, and x and probability. | |
223 | // Uses for these parameter estimators are unclear. | |
224 | ||
225 | static RealType find_alpha( | |
226 | RealType beta, // from beta. | |
227 | RealType x, // x. | |
228 | RealType probability) // cdf | |
229 | { | |
230 | static const char* function = "boost::math::beta_distribution<%1%>::find_alpha"; | |
231 | RealType result = 0; // of error checks. | |
232 | if(false == | |
233 | ( | |
234 | beta_detail::check_prob(function, probability, &result, Policy()) | |
235 | && | |
236 | beta_detail::check_beta(function, beta, &result, Policy()) | |
237 | && | |
238 | beta_detail::check_x(function, x, &result, Policy()) | |
239 | ) | |
240 | ) | |
241 | { | |
242 | return result; | |
243 | } | |
244 | return ibeta_inva(beta, x, probability, Policy()); | |
245 | } // RealType find_alpha(beta, a, probability) | |
246 | ||
247 | static RealType find_beta( | |
248 | // ibeta_invb(T b, T x, T p); (alpha, x, cdf,) | |
249 | RealType alpha, // alpha. | |
250 | RealType x, // probability x. | |
251 | RealType probability) // probability cdf. | |
252 | { | |
253 | static const char* function = "boost::math::beta_distribution<%1%>::find_beta"; | |
254 | RealType result = 0; // of error checks. | |
255 | if(false == | |
256 | ( | |
257 | beta_detail::check_prob(function, probability, &result, Policy()) | |
258 | && | |
259 | beta_detail::check_alpha(function, alpha, &result, Policy()) | |
260 | && | |
261 | beta_detail::check_x(function, x, &result, Policy()) | |
262 | ) | |
263 | ) | |
264 | { | |
265 | return result; | |
266 | } | |
267 | return ibeta_invb(alpha, x, probability, Policy()); | |
268 | } // RealType find_beta(alpha, x, probability) | |
269 | ||
270 | private: | |
271 | RealType m_alpha; // Two parameters of the beta distribution. | |
272 | RealType m_beta; | |
273 | }; // template <class RealType, class Policy> class beta_distribution | |
274 | ||
275 | template <class RealType, class Policy> | |
276 | inline const std::pair<RealType, RealType> range(const beta_distribution<RealType, Policy>& /* dist */) | |
277 | { // Range of permissible values for random variable x. | |
278 | using boost::math::tools::max_value; | |
279 | return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1)); | |
280 | } | |
281 | ||
282 | template <class RealType, class Policy> | |
283 | inline const std::pair<RealType, RealType> support(const beta_distribution<RealType, Policy>& /* dist */) | |
284 | { // Range of supported values for random variable x. | |
285 | // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. | |
286 | return std::pair<RealType, RealType>(static_cast<RealType>(0), static_cast<RealType>(1)); | |
287 | } | |
288 | ||
289 | template <class RealType, class Policy> | |
290 | inline RealType mean(const beta_distribution<RealType, Policy>& dist) | |
291 | { // Mean of beta distribution = np. | |
292 | return dist.alpha() / (dist.alpha() + dist.beta()); | |
293 | } // mean | |
294 | ||
295 | template <class RealType, class Policy> | |
296 | inline RealType variance(const beta_distribution<RealType, Policy>& dist) | |
297 | { // Variance of beta distribution = np(1-p). | |
298 | RealType a = dist.alpha(); | |
299 | RealType b = dist.beta(); | |
300 | return (a * b) / ((a + b ) * (a + b) * (a + b + 1)); | |
301 | } // variance | |
302 | ||
303 | template <class RealType, class Policy> | |
304 | inline RealType mode(const beta_distribution<RealType, Policy>& dist) | |
305 | { | |
306 | static const char* function = "boost::math::mode(beta_distribution<%1%> const&)"; | |
307 | ||
308 | RealType result; | |
309 | if ((dist.alpha() <= 1)) | |
310 | { | |
311 | result = policies::raise_domain_error<RealType>( | |
312 | function, | |
313 | "mode undefined for alpha = %1%, must be > 1!", dist.alpha(), Policy()); | |
314 | return result; | |
315 | } | |
316 | ||
317 | if ((dist.beta() <= 1)) | |
318 | { | |
319 | result = policies::raise_domain_error<RealType>( | |
320 | function, | |
321 | "mode undefined for beta = %1%, must be > 1!", dist.beta(), Policy()); | |
322 | return result; | |
323 | } | |
324 | RealType a = dist.alpha(); | |
325 | RealType b = dist.beta(); | |
326 | return (a-1) / (a + b - 2); | |
327 | } // mode | |
328 | ||
329 | //template <class RealType, class Policy> | |
330 | //inline RealType median(const beta_distribution<RealType, Policy>& dist) | |
331 | //{ // Median of beta distribution is not defined. | |
332 | // return tools::domain_error<RealType>(function, "Median is not implemented, result is %1%!", std::numeric_limits<RealType>::quiet_NaN()); | |
333 | //} // median | |
334 | ||
335 | //But WILL be provided by the derived accessor as quantile(0.5). | |
336 | ||
337 | template <class RealType, class Policy> | |
338 | inline RealType skewness(const beta_distribution<RealType, Policy>& dist) | |
339 | { | |
340 | BOOST_MATH_STD_USING // ADL of std functions. | |
341 | RealType a = dist.alpha(); | |
342 | RealType b = dist.beta(); | |
343 | return (2 * (b-a) * sqrt(a + b + 1)) / ((a + b + 2) * sqrt(a * b)); | |
344 | } // skewness | |
345 | ||
346 | template <class RealType, class Policy> | |
347 | inline RealType kurtosis_excess(const beta_distribution<RealType, Policy>& dist) | |
348 | { | |
349 | RealType a = dist.alpha(); | |
350 | RealType b = dist.beta(); | |
351 | RealType a_2 = a * a; | |
352 | RealType n = 6 * (a_2 * a - a_2 * (2 * b - 1) + b * b * (b + 1) - 2 * a * b * (b + 2)); | |
353 | RealType d = a * b * (a + b + 2) * (a + b + 3); | |
354 | return n / d; | |
355 | } // kurtosis_excess | |
356 | ||
357 | template <class RealType, class Policy> | |
358 | inline RealType kurtosis(const beta_distribution<RealType, Policy>& dist) | |
359 | { | |
360 | return 3 + kurtosis_excess(dist); | |
361 | } // kurtosis | |
362 | ||
363 | template <class RealType, class Policy> | |
364 | inline RealType pdf(const beta_distribution<RealType, Policy>& dist, const RealType& x) | |
365 | { // Probability Density/Mass Function. | |
366 | BOOST_FPU_EXCEPTION_GUARD | |
367 | ||
368 | static const char* function = "boost::math::pdf(beta_distribution<%1%> const&, %1%)"; | |
369 | ||
370 | BOOST_MATH_STD_USING // for ADL of std functions | |
371 | ||
372 | RealType a = dist.alpha(); | |
373 | RealType b = dist.beta(); | |
374 | ||
375 | // Argument checks: | |
376 | RealType result = 0; | |
377 | if(false == beta_detail::check_dist_and_x( | |
378 | function, | |
379 | a, b, x, | |
380 | &result, Policy())) | |
381 | { | |
382 | return result; | |
383 | } | |
384 | using boost::math::beta; | |
385 | return ibeta_derivative(a, b, x, Policy()); | |
386 | ||
387 | ||
388 | template <class RealType, class Policy> | |
389 | inline RealType cdf(const beta_distribution<RealType, Policy>& dist, const RealType& x) | |
390 | { // Cumulative Distribution Function beta. | |
391 | BOOST_MATH_STD_USING // for ADL of std functions | |
392 | ||
393 | static const char* function = "boost::math::cdf(beta_distribution<%1%> const&, %1%)"; | |
394 | ||
395 | RealType a = dist.alpha(); | |
396 | RealType b = dist.beta(); | |
397 | ||
398 | // Argument checks: | |
399 | RealType result = 0; | |
400 | if(false == beta_detail::check_dist_and_x( | |
401 | function, | |
402 | a, b, x, | |
403 | &result, Policy())) | |
404 | { | |
405 | return result; | |
406 | } | |
407 | // Special cases: | |
408 | if (x == 0) | |
409 | { | |
410 | return 0; | |
411 | } | |
412 | else if (x == 1) | |
413 | { | |
414 | return 1; | |
415 | } | |
416 | return ibeta(a, b, x, Policy()); | |
417 | } // beta cdf | |
418 | ||
419 | template <class RealType, class Policy> | |
420 | inline RealType cdf(const complemented2_type<beta_distribution<RealType, Policy>, RealType>& c) | |
421 | { // Complemented Cumulative Distribution Function beta. | |
422 | ||
423 | BOOST_MATH_STD_USING // for ADL of std functions | |
424 | ||
425 | static const char* function = "boost::math::cdf(beta_distribution<%1%> const&, %1%)"; | |
426 | ||
427 | RealType const& x = c.param; | |
428 | beta_distribution<RealType, Policy> const& dist = c.dist; | |
429 | RealType a = dist.alpha(); | |
430 | RealType b = dist.beta(); | |
431 | ||
432 | // Argument checks: | |
433 | RealType result = 0; | |
434 | if(false == beta_detail::check_dist_and_x( | |
435 | function, | |
436 | a, b, x, | |
437 | &result, Policy())) | |
438 | { | |
439 | return result; | |
440 | } | |
441 | if (x == 0) | |
442 | { | |
443 | return 1; | |
444 | } | |
445 | else if (x == 1) | |
446 | { | |
447 | return 0; | |
448 | } | |
449 | // Calculate cdf beta using the incomplete beta function. | |
450 | // Use of ibeta here prevents cancellation errors in calculating | |
451 | // 1 - x if x is very small, perhaps smaller than machine epsilon. | |
452 | return ibetac(a, b, x, Policy()); | |
453 | } // beta cdf | |
454 | ||
455 | template <class RealType, class Policy> | |
456 | inline RealType quantile(const beta_distribution<RealType, Policy>& dist, const RealType& p) | |
457 | { // Quantile or Percent Point beta function or | |
458 | // Inverse Cumulative probability distribution function CDF. | |
459 | // Return x (0 <= x <= 1), | |
460 | // for a given probability p (0 <= p <= 1). | |
461 | // These functions take a probability as an argument | |
462 | // and return a value such that the probability that a random variable x | |
463 | // will be less than or equal to that value | |
464 | // is whatever probability you supplied as an argument. | |
465 | ||
466 | static const char* function = "boost::math::quantile(beta_distribution<%1%> const&, %1%)"; | |
467 | ||
468 | RealType result = 0; // of argument checks: | |
469 | RealType a = dist.alpha(); | |
470 | RealType b = dist.beta(); | |
471 | if(false == beta_detail::check_dist_and_prob( | |
472 | function, | |
473 | a, b, p, | |
474 | &result, Policy())) | |
475 | { | |
476 | return result; | |
477 | } | |
478 | // Special cases: | |
479 | if (p == 0) | |
480 | { | |
481 | return 0; | |
482 | } | |
483 | if (p == 1) | |
484 | { | |
485 | return 1; | |
486 | } | |
487 | return ibeta_inv(a, b, p, static_cast<RealType*>(0), Policy()); | |
488 | } // quantile | |
489 | ||
490 | template <class RealType, class Policy> | |
491 | inline RealType quantile(const complemented2_type<beta_distribution<RealType, Policy>, RealType>& c) | |
492 | { // Complement Quantile or Percent Point beta function . | |
493 | // Return the number of expected x for a given | |
494 | // complement of the probability q. | |
495 | ||
496 | static const char* function = "boost::math::quantile(beta_distribution<%1%> const&, %1%)"; | |
497 | ||
498 | // | |
499 | // Error checks: | |
500 | RealType q = c.param; | |
501 | const beta_distribution<RealType, Policy>& dist = c.dist; | |
502 | RealType result = 0; | |
503 | RealType a = dist.alpha(); | |
504 | RealType b = dist.beta(); | |
505 | if(false == beta_detail::check_dist_and_prob( | |
506 | function, | |
507 | a, | |
508 | b, | |
509 | q, | |
510 | &result, Policy())) | |
511 | { | |
512 | return result; | |
513 | } | |
514 | // Special cases: | |
515 | if(q == 1) | |
516 | { | |
517 | return 0; | |
518 | } | |
519 | if(q == 0) | |
520 | { | |
521 | return 1; | |
522 | } | |
523 | ||
524 | return ibetac_inv(a, b, q, static_cast<RealType*>(0), Policy()); | |
525 | } // Quantile Complement | |
526 | ||
527 | } // namespace math | |
528 | } // namespace boost | |
529 | ||
530 | // This include must be at the end, *after* the accessors | |
531 | // for this distribution have been defined, in order to | |
532 | // keep compilers that support two-phase lookup happy. | |
533 | #include <boost/math/distributions/detail/derived_accessors.hpp> | |
534 | ||
535 | #if defined (BOOST_MSVC) | |
536 | # pragma warning(pop) | |
537 | #endif | |
538 | ||
539 | #endif // BOOST_MATH_DIST_BETA_HPP | |
540 | ||
541 |