1// boost\math\distributions\binomial.hpp
2
3// Copyright John Maddock 2006.
4// Copyright Paul A. Bristow 2007.
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/binomial_distribution
12
13// Binomial distribution is the discrete probability distribution of
14// the number (k) of successes, in a sequence of
15// n independent (yes or no, success or failure) Bernoulli trials.
16
17// It expresses the probability of a number of events occurring in a fixed time
18// if these events occur with a known average rate (probability of success),
19// and are independent of the time since the last event.
20
21// The number of cars that pass through a certain point on a road during a given period of time.
22// The number of spelling mistakes a secretary makes while typing a single page.
23// The number of phone calls at a call center per minute.
24// The number of times a web server is accessed per minute.
25// The number of light bulbs that burn out in a certain amount of time.
26// The number of roadkill found per unit length of road
27
28// http://en.wikipedia.org/wiki/binomial_distribution
29
30// Given a sample of N measured values k[i],
31// we wish to estimate the value of the parameter x (mean)
32// of the binomial population from which the sample was drawn.
33// To calculate the maximum likelihood value = 1/N sum i = 1 to N of k[i]
34
35// Also may want a function for EXACTLY k.
36
37// And probability that there are EXACTLY k occurrences is
38// exp(-x) * pow(x, k) / factorial(k)
39// where x is expected occurrences (mean) during the given interval.
40// For example, if events occur, on average, every 4 min,
41// and we are interested in number of events occurring in 10 min,
42// then x = 10/4 = 2.5
43
44// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366i.htm
45
46// The binomial distribution is used when there are
47// exactly two mutually exclusive outcomes of a trial.
48// These outcomes are appropriately labeled "success" and "failure".
49// The binomial distribution is used to obtain
50// the probability of observing x successes in N trials,
51// with the probability of success on a single trial denoted by p.
52// The binomial distribution assumes that p is fixed for all trials.
53
54// P(x, p, n) = n!/(x! * (n-x)!) * p^x * (1-p)^(n-x)
55
56// http://mathworld.wolfram.com/BinomialCoefficient.html
57
58// The binomial coefficient (n; k) is the number of ways of picking
59// k unordered outcomes from n possibilities,
60// also known as a combination or combinatorial number.
61// The symbols _nC_k and (n; k) are used to denote a binomial coefficient,
62// and are sometimes read as "n choose k."
63// (n; k) therefore gives the number of k-subsets possible out of a set of n distinct items.
64
65// For example:
66// The 2-subsets of {1,2,3,4} are the six pairs {1,2}, {1,3}, {1,4}, {2,3}, {2,4}, and {3,4}, so (4; 2)==6.
67
68// http://functions.wolfram.com/GammaBetaErf/Binomial/ for evaluation.
69
70// But note that the binomial distribution
71// (like others including the poisson, negative binomial & Bernoulli)
72// is strictly defined as a discrete function: only integral values of k are envisaged.
73// However because of the method of calculation using a continuous gamma function,
74// it is convenient to treat it as if a continuous function,
75// and permit non-integral values of k.
76// To enforce the strict mathematical model, users should use floor or ceil functions
77// on k outside this function to ensure that k is integral.
78
79#ifndef BOOST_MATH_SPECIAL_BINOMIAL_HPP
80#define BOOST_MATH_SPECIAL_BINOMIAL_HPP
81
82#include <boost/math/distributions/fwd.hpp>
83#include <boost/math/special_functions/beta.hpp> // for incomplete beta.
84#include <boost/math/distributions/complement.hpp> // complements
85#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
86#include <boost/math/distributions/detail/inv_discrete_quantile.hpp> // error checks
87#include <boost/math/special_functions/fpclassify.hpp> // isnan.
88#include <boost/math/tools/roots.hpp> // for root finding.
89
90#include <utility>
91
92namespace boost
93{
94 namespace math
95 {
96
97 template <class RealType, class Policy>
98 class binomial_distribution;
99
100 namespace binomial_detail{
101 // common error checking routines for binomial distribution functions:
102 template <class RealType, class Policy>
103 inline bool check_N(const char* function, const RealType& N, RealType* result, const Policy& pol)
104 {
105 if((N < 0) || !(boost::math::isfinite)(N))
106 {
107 *result = policies::raise_domain_error<RealType>(
108 function,
109 "Number of Trials argument is %1%, but must be >= 0 !", N, pol);
110 return false;
111 }
112 return true;
113 }
114 template <class RealType, class Policy>
115 inline bool check_success_fraction(const char* function, const RealType& p, RealType* result, const Policy& pol)
116 {
117 if((p < 0) || (p > 1) || !(boost::math::isfinite)(p))
118 {
119 *result = policies::raise_domain_error<RealType>(
120 function,
121 "Success fraction argument is %1%, but must be >= 0 and <= 1 !", p, pol);
122 return false;
123 }
124 return true;
125 }
126 template <class RealType, class Policy>
127 inline bool check_dist(const char* function, const RealType& N, const RealType& p, RealType* result, const Policy& pol)
128 {
129 return check_success_fraction(
130 function, p, result, pol)
131 && check_N(
132 function, N, result, pol);
133 }
134 template <class RealType, class Policy>
135 inline bool check_dist_and_k(const char* function, const RealType& N, const RealType& p, RealType k, RealType* result, const Policy& pol)
136 {
137 if(check_dist(function, N, p, result, pol) == false)
138 return false;
139 if((k < 0) || !(boost::math::isfinite)(k))
140 {
141 *result = policies::raise_domain_error<RealType>(
142 function,
143 "Number of Successes argument is %1%, but must be >= 0 !", k, pol);
144 return false;
145 }
146 if(k > N)
147 {
148 *result = policies::raise_domain_error<RealType>(
149 function,
150 "Number of Successes argument is %1%, but must be <= Number of Trials !", k, pol);
151 return false;
152 }
153 return true;
154 }
155 template <class RealType, class Policy>
156 inline bool check_dist_and_prob(const char* function, const RealType& N, RealType p, RealType prob, RealType* result, const Policy& pol)
157 {
158 if((check_dist(function, N, p, result, pol) && detail::check_probability(function, prob, result, pol)) == false)
159 return false;
160 return true;
161 }
162
163 template <class T, class Policy>
164 T inverse_binomial_cornish_fisher(T n, T sf, T p, T q, const Policy& pol)
165 {
166 BOOST_MATH_STD_USING
167 // mean:
168 T m = n * sf;
169 // standard deviation:
170 T sigma = sqrt(n * sf * (1 - sf));
171 // skewness
172 T sk = (1 - 2 * sf) / sigma;
173 // kurtosis:
174 // T k = (1 - 6 * sf * (1 - sf) ) / (n * sf * (1 - sf));
175 // Get the inverse of a std normal distribution:
176 T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();
177 // Set the sign:
178 if(p < 0.5)
179 x = -x;
180 T x2 = x * x;
181 // w is correction term due to skewness
182 T w = x + sk * (x2 - 1) / 6;
183 /*
184 // Add on correction due to kurtosis.
185 // Disabled for now, seems to make things worse?
186 //
187 if(n >= 10)
188 w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
189 */
190 w = m + sigma * w;
191 if(w < tools::min_value<T>())
192 return sqrt(tools::min_value<T>());
193 if(w > n)
194 return n;
195 return w;
196 }
197
198 template <class RealType, class Policy>
199 RealType quantile_imp(const binomial_distribution<RealType, Policy>& dist, const RealType& p, const RealType& q, bool comp)
200 { // Quantile or Percent Point Binomial function.
201 // Return the number of expected successes k,
202 // for a given probability p.
203 //
204 // Error checks:
205 BOOST_MATH_STD_USING // ADL of std names
206 RealType result = 0;
207 RealType trials = dist.trials();
208 RealType success_fraction = dist.success_fraction();
209 if(false == binomial_detail::check_dist_and_prob(
210 "boost::math::quantile(binomial_distribution<%1%> const&, %1%)",
211 trials,
212 success_fraction,
213 p,
214 &result, Policy()))
215 {
216 return result;
217 }
218
219 // Special cases:
220 //
221 if(p == 0)
222 { // There may actually be no answer to this question,
223 // since the probability of zero successes may be non-zero,
224 // but zero is the best we can do:
225 return 0;
226 }
227 if(p == 1 || success_fraction == 1)
228 { // Probability of n or fewer successes is always one,
229 // so n is the most sensible answer here:
230 return trials;
231 }
232 if (p <= pow(1 - success_fraction, trials))
233 { // p <= pdf(dist, 0) == cdf(dist, 0)
234 return 0; // So the only reasonable result is zero.
235 } // And root finder would fail otherwise.
236
237 // Solve for quantile numerically:
238 //
239 RealType guess = binomial_detail::inverse_binomial_cornish_fisher(trials, success_fraction, p, q, Policy());
240 RealType factor = 8;
241 if(trials > 100)
242 factor = 1.01f; // guess is pretty accurate
243 else if((trials > 10) && (trials - 1 > guess) && (guess > 3))
244 factor = 1.15f; // less accurate but OK.
245 else if(trials < 10)
246 {
247 // pretty inaccurate guess in this area:
248 if(guess > trials / 64)
249 {
250 guess = trials / 4;
251 factor = 2;
252 }
253 else
254 guess = trials / 1024;
255 }
256 else
257 factor = 2; // trials largish, but in far tails.
258
259 typedef typename Policy::discrete_quantile_type discrete_quantile_type;
260 std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
261 result = detail::inverse_discrete_quantile(
262 dist,
263 comp ? q : p,
264 comp,
265 guess,
266 factor,
267 RealType(1),
268 discrete_quantile_type(),
269 max_iter);
270 return result;
271 } // quantile
272
273 }
274
275 template <class RealType = double, class Policy = policies::policy<> >
276 class binomial_distribution
277 {
278 public:
279 typedef RealType value_type;
280 typedef Policy policy_type;
281
282 binomial_distribution(RealType n = 1, RealType p = 0.5) : m_n(n), m_p(p)
283 { // Default n = 1 is the Bernoulli distribution
284 // with equal probability of 'heads' or 'tails.
285 RealType r;
286 binomial_detail::check_dist(
287 "boost::math::binomial_distribution<%1%>::binomial_distribution",
288 m_n,
289 m_p,
290 &r, Policy());
291 } // binomial_distribution constructor.
292
293 RealType success_fraction() const
294 { // Probability.
295 return m_p;
296 }
297 RealType trials() const
298 { // Total number of trials.
299 return m_n;
300 }
301
302 enum interval_type{
303 clopper_pearson_exact_interval,
304 jeffreys_prior_interval
305 };
306
307 //
308 // Estimation of the success fraction parameter.
309 // The best estimate is actually simply successes/trials,
310 // these functions are used
311 // to obtain confidence intervals for the success fraction.
312 //
313 static RealType find_lower_bound_on_p(
314 RealType trials,
315 RealType successes,
316 RealType probability,
317 interval_type t = clopper_pearson_exact_interval)
318 {
319 static const char* function = "boost::math::binomial_distribution<%1%>::find_lower_bound_on_p";
320 // Error checks:
321 RealType result = 0;
322 if(false == binomial_detail::check_dist_and_k(
323 function, trials, RealType(0), successes, &result, Policy())
324 &&
325 binomial_detail::check_dist_and_prob(
326 function, trials, RealType(0), probability, &result, Policy()))
327 { return result; }
328
329 if(successes == 0)
330 return 0;
331
332 // NOTE!!! The Clopper Pearson formula uses "successes" not
333 // "successes+1" as usual to get the lower bound,
334 // see http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm
335 return (t == clopper_pearson_exact_interval) ? ibeta_inv(successes, trials - successes + 1, probability, static_cast<RealType*>(nullptr), Policy())
336 : ibeta_inv(successes + 0.5f, trials - successes + 0.5f, probability, static_cast<RealType*>(nullptr), Policy());
337 }
338 static RealType find_upper_bound_on_p(
339 RealType trials,
340 RealType successes,
341 RealType probability,
342 interval_type t = clopper_pearson_exact_interval)
343 {
344 static const char* function = "boost::math::binomial_distribution<%1%>::find_upper_bound_on_p";
345 // Error checks:
346 RealType result = 0;
347 if(false == binomial_detail::check_dist_and_k(
348 function, trials, RealType(0), successes, &result, Policy())
349 &&
350 binomial_detail::check_dist_and_prob(
351 function, trials, RealType(0), probability, &result, Policy()))
352 { return result; }
353
354 if(trials == successes)
355 return 1;
356
357 return (t == clopper_pearson_exact_interval) ? ibetac_inv(successes + 1, trials - successes, probability, static_cast<RealType*>(nullptr), Policy())
358 : ibetac_inv(successes + 0.5f, trials - successes + 0.5f, probability, static_cast<RealType*>(nullptr), Policy());
359 }
360 // Estimate number of trials parameter:
361 //
362 // "How many trials do I need to be P% sure of seeing k events?"
363 // or
364 // "How many trials can I have to be P% sure of seeing fewer than k events?"
365 //
366 static RealType find_minimum_number_of_trials(
367 RealType k, // number of events
368 RealType p, // success fraction
369 RealType alpha) // risk level
370 {
371 static const char* function = "boost::math::binomial_distribution<%1%>::find_minimum_number_of_trials";
372 // Error checks:
373 RealType result = 0;
374 if(false == binomial_detail::check_dist_and_k(
375 function, k, p, k, &result, Policy())
376 &&
377 binomial_detail::check_dist_and_prob(
378 function, k, p, alpha, &result, Policy()))
379 { return result; }
380
381 result = ibetac_invb(k + 1, p, alpha, Policy()); // returns n - k
382 return result + k;
383 }
384
385 static RealType find_maximum_number_of_trials(
386 RealType k, // number of events
387 RealType p, // success fraction
388 RealType alpha) // risk level
389 {
390 static const char* function = "boost::math::binomial_distribution<%1%>::find_maximum_number_of_trials";
391 // Error checks:
392 RealType result = 0;
393 if(false == binomial_detail::check_dist_and_k(
394 function, k, p, k, &result, Policy())
395 &&
396 binomial_detail::check_dist_and_prob(
397 function, k, p, alpha, &result, Policy()))
398 { return result; }
399
400 result = ibeta_invb(k + 1, p, alpha, Policy()); // returns n - k
401 return result + k;
402 }
403
404 private:
405 RealType m_n; // Not sure if this shouldn't be an int?
406 RealType m_p; // success_fraction
407 }; // template <class RealType, class Policy> class binomial_distribution
408
409 typedef binomial_distribution<> binomial;
410 // typedef binomial_distribution<double> binomial;
411 // IS now included since no longer a name clash with function binomial.
412 //typedef binomial_distribution<double> binomial; // Reserved name of type double.
413
414 #ifdef __cpp_deduction_guides
415 template <class RealType>
416 binomial_distribution(RealType)->binomial_distribution<typename boost::math::tools::promote_args<RealType>::type>;
417 template <class RealType>
418 binomial_distribution(RealType,RealType)->binomial_distribution<typename boost::math::tools::promote_args<RealType>::type>;
419 #endif
420
421 template <class RealType, class Policy>
422 const std::pair<RealType, RealType> range(const binomial_distribution<RealType, Policy>& dist)
423 { // Range of permissible values for random variable k.
424 using boost::math::tools::max_value;
425 return std::pair<RealType, RealType>(static_cast<RealType>(0), dist.trials());
426 }
427
428 template <class RealType, class Policy>
429 const std::pair<RealType, RealType> support(const binomial_distribution<RealType, Policy>& dist)
430 { // Range of supported values for random variable k.
431 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
432 return std::pair<RealType, RealType>(static_cast<RealType>(0), dist.trials());
433 }
434
435 template <class RealType, class Policy>
436 inline RealType mean(const binomial_distribution<RealType, Policy>& dist)
437 { // Mean of Binomial distribution = np.
438 return dist.trials() * dist.success_fraction();
439 } // mean
440
441 template <class RealType, class Policy>
442 inline RealType variance(const binomial_distribution<RealType, Policy>& dist)
443 { // Variance of Binomial distribution = np(1-p).
444 return dist.trials() * dist.success_fraction() * (1 - dist.success_fraction());
445 } // variance
446
447 template <class RealType, class Policy>
448 RealType pdf(const binomial_distribution<RealType, Policy>& dist, const RealType& k)
449 { // Probability Density/Mass Function.
450 BOOST_FPU_EXCEPTION_GUARD
451
452 BOOST_MATH_STD_USING // for ADL of std functions
453
454 RealType n = dist.trials();
455
456 // Error check:
457 RealType result = 0; // initialization silences some compiler warnings
458 if(false == binomial_detail::check_dist_and_k(
459 "boost::math::pdf(binomial_distribution<%1%> const&, %1%)",
460 n,
461 dist.success_fraction(),
462 k,
463 &result, Policy()))
464 {
465 return result;
466 }
467
468 // Special cases of success_fraction, regardless of k successes and regardless of n trials.
469 if (dist.success_fraction() == 0)
470 { // probability of zero successes is 1:
471 return static_cast<RealType>(k == 0 ? 1 : 0);
472 }
473 if (dist.success_fraction() == 1)
474 { // probability of n successes is 1:
475 return static_cast<RealType>(k == n ? 1 : 0);
476 }
477 // k argument may be integral, signed, or unsigned, or floating point.
478 // If necessary, it has already been promoted from an integral type.
479 if (n == 0)
480 {
481 return 1; // Probability = 1 = certainty.
482 }
483 if (k == n)
484 { // binomial coeffic (n n) = 1,
485 // n ^ 0 = 1
486 return pow(dist.success_fraction(), k); // * pow((1 - dist.success_fraction()), (n - k)) = 1
487 }
488
489 // Probability of getting exactly k successes
490 // if C(n, k) is the binomial coefficient then:
491 //
492 // f(k; n,p) = C(n, k) * p^k * (1-p)^(n-k)
493 // = (n!/(k!(n-k)!)) * p^k * (1-p)^(n-k)
494 // = (tgamma(n+1) / (tgamma(k+1)*tgamma(n-k+1))) * p^k * (1-p)^(n-k)
495 // = p^k (1-p)^(n-k) / (beta(k+1, n-k+1) * (n+1))
496 // = ibeta_derivative(k+1, n-k+1, p) / (n+1)
497 //
498 using boost::math::ibeta_derivative; // a, b, x
499 return ibeta_derivative(k+1, n-k+1, dist.success_fraction(), Policy()) / (n+1);
500
501 } // pdf
502
503 template <class RealType, class Policy>
504 inline RealType cdf(const binomial_distribution<RealType, Policy>& dist, const RealType& k)
505 { // Cumulative Distribution Function Binomial.
506 // The random variate k is the number of successes in n trials.
507 // k argument may be integral, signed, or unsigned, or floating point.
508 // If necessary, it has already been promoted from an integral type.
509
510 // Returns the sum of the terms 0 through k of the Binomial Probability Density/Mass:
511 //
512 // i=k
513 // -- ( n ) i n-i
514 // > | | p (1-p)
515 // -- ( i )
516 // i=0
517
518 // The terms are not summed directly instead
519 // the incomplete beta integral is employed,
520 // according to the formula:
521 // P = I[1-p]( n-k, k+1).
522 // = 1 - I[p](k + 1, n - k)
523
524 BOOST_MATH_STD_USING // for ADL of std functions
525
526 RealType n = dist.trials();
527 RealType p = dist.success_fraction();
528
529 // Error check:
530 RealType result = 0;
531 if(false == binomial_detail::check_dist_and_k(
532 "boost::math::cdf(binomial_distribution<%1%> const&, %1%)",
533 n,
534 p,
535 k,
536 &result, Policy()))
537 {
538 return result;
539 }
540 if (k == n)
541 {
542 return 1;
543 }
544
545 // Special cases, regardless of k.
546 if (p == 0)
547 { // This need explanation:
548 // the pdf is zero for all cases except when k == 0.
549 // For zero p the probability of zero successes is one.
550 // Therefore the cdf is always 1:
551 // the probability of k or *fewer* successes is always 1
552 // if there are never any successes!
553 return 1;
554 }
555 if (p == 1)
556 { // This is correct but needs explanation:
557 // when k = 1
558 // all the cdf and pdf values are zero *except* when k == n,
559 // and that case has been handled above already.
560 return 0;
561 }
562 //
563 // P = I[1-p](n - k, k + 1)
564 // = 1 - I[p](k + 1, n - k)
565 // Use of ibetac here prevents cancellation errors in calculating
566 // 1-p if p is very small, perhaps smaller than machine epsilon.
567 //
568 // Note that we do not use a finite sum here, since the incomplete
569 // beta uses a finite sum internally for integer arguments, so
570 // we'll just let it take care of the necessary logic.
571 //
572 return ibetac(k + 1, n - k, p, Policy());
573 } // binomial cdf
574
575 template <class RealType, class Policy>
576 inline RealType cdf(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c)
577 { // Complemented Cumulative Distribution Function Binomial.
578 // The random variate k is the number of successes in n trials.
579 // k argument may be integral, signed, or unsigned, or floating point.
580 // If necessary, it has already been promoted from an integral type.
581
582 // Returns the sum of the terms k+1 through n of the Binomial Probability Density/Mass:
583 //
584 // i=n
585 // -- ( n ) i n-i
586 // > | | p (1-p)
587 // -- ( i )
588 // i=k+1
589
590 // The terms are not summed directly instead
591 // the incomplete beta integral is employed,
592 // according to the formula:
593 // Q = 1 -I[1-p]( n-k, k+1).
594 // = I[p](k + 1, n - k)
595
596 BOOST_MATH_STD_USING // for ADL of std functions
597
598 RealType const& k = c.param;
599 binomial_distribution<RealType, Policy> const& dist = c.dist;
600 RealType n = dist.trials();
601 RealType p = dist.success_fraction();
602
603 // Error checks:
604 RealType result = 0;
605 if(false == binomial_detail::check_dist_and_k(
606 "boost::math::cdf(binomial_distribution<%1%> const&, %1%)",
607 n,
608 p,
609 k,
610 &result, Policy()))
611 {
612 return result;
613 }
614
615 if (k == n)
616 { // Probability of greater than n successes is necessarily zero:
617 return 0;
618 }
619
620 // Special cases, regardless of k.
621 if (p == 0)
622 {
623 // This need explanation: the pdf is zero for all
624 // cases except when k == 0. For zero p the probability
625 // of zero successes is one. Therefore the cdf is always
626 // 1: the probability of *more than* k successes is always 0
627 // if there are never any successes!
628 return 0;
629 }
630 if (p == 1)
631 {
632 // This needs explanation, when p = 1
633 // we always have n successes, so the probability
634 // of more than k successes is 1 as long as k < n.
635 // The k == n case has already been handled above.
636 return 1;
637 }
638 //
639 // Calculate cdf binomial using the incomplete beta function.
640 // Q = 1 -I[1-p](n - k, k + 1)
641 // = I[p](k + 1, n - k)
642 // Use of ibeta here prevents cancellation errors in calculating
643 // 1-p if p is very small, perhaps smaller than machine epsilon.
644 //
645 // Note that we do not use a finite sum here, since the incomplete
646 // beta uses a finite sum internally for integer arguments, so
647 // we'll just let it take care of the necessary logic.
648 //
649 return ibeta(k + 1, n - k, p, Policy());
650 } // binomial cdf
651
652 template <class RealType, class Policy>
653 inline RealType quantile(const binomial_distribution<RealType, Policy>& dist, const RealType& p)
654 {
655 return binomial_detail::quantile_imp(dist, p, RealType(1-p), false);
656 } // quantile
657
658 template <class RealType, class Policy>
659 RealType quantile(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c)
660 {
661 return binomial_detail::quantile_imp(c.dist, RealType(1-c.param), c.param, true);
662 } // quantile
663
664 template <class RealType, class Policy>
665 inline RealType mode(const binomial_distribution<RealType, Policy>& dist)
666 {
667 BOOST_MATH_STD_USING // ADL of std functions.
668 RealType p = dist.success_fraction();
669 RealType n = dist.trials();
670 return floor(p * (n + 1));
671 }
672
673 template <class RealType, class Policy>
674 inline RealType median(const binomial_distribution<RealType, Policy>& dist)
675 { // Bounds for the median of the negative binomial distribution
676 // VAN DE VEN R. ; WEBER N. C. ;
677 // Univ. Sydney, school mathematics statistics, Sydney N.S.W. 2006, AUSTRALIE
678 // Metrika (Metrika) ISSN 0026-1335 CODEN MTRKA8
679 // 1993, vol. 40, no3-4, pp. 185-189 (4 ref.)
680
681 // Bounds for median and 50 percentage point of binomial and negative binomial distribution
682 // Metrika, ISSN 0026-1335 (Print) 1435-926X (Online)
683 // Volume 41, Number 1 / December, 1994, DOI 10.1007/BF01895303
684 BOOST_MATH_STD_USING // ADL of std functions.
685 RealType p = dist.success_fraction();
686 RealType n = dist.trials();
687 // Wikipedia says one of floor(np) -1, floor (np), floor(np) +1
688 return floor(p * n); // Chose the middle value.
689 }
690
691 template <class RealType, class Policy>
692 inline RealType skewness(const binomial_distribution<RealType, Policy>& dist)
693 {
694 BOOST_MATH_STD_USING // ADL of std functions.
695 RealType p = dist.success_fraction();
696 RealType n = dist.trials();
697 return (1 - 2 * p) / sqrt(n * p * (1 - p));
698 }
699
700 template <class RealType, class Policy>
701 inline RealType kurtosis(const binomial_distribution<RealType, Policy>& dist)
702 {
703 RealType p = dist.success_fraction();
704 RealType n = dist.trials();
705 return 3 - 6 / n + 1 / (n * p * (1 - p));
706 }
707
708 template <class RealType, class Policy>
709 inline RealType kurtosis_excess(const binomial_distribution<RealType, Policy>& dist)
710 {
711 RealType p = dist.success_fraction();
712 RealType q = 1 - p;
713 RealType n = dist.trials();
714 return (1 - 6 * p * q) / (n * p * q);
715 }
716
717 } // namespace math
718 } // namespace boost
719
720// This include must be at the end, *after* the accessors
721// for this distribution have been defined, in order to
722// keep compilers that support two-phase lookup happy.
723#include <boost/math/distributions/detail/derived_accessors.hpp>
724
725#endif // BOOST_MATH_SPECIAL_BINOMIAL_HPP
726
727
728

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