]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 continous 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 | ||
92 | namespace 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) | |
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 | if(success_fraction == 1) | |
237 | { // our formulae break down in this case: | |
238 | return p > 0.5f ? trials : 0; | |
239 | } | |
240 | ||
241 | // Solve for quantile numerically: | |
242 | // | |
243 | RealType guess = binomial_detail::inverse_binomial_cornish_fisher(trials, success_fraction, p, q, Policy()); | |
244 | RealType factor = 8; | |
245 | if(trials > 100) | |
246 | factor = 1.01f; // guess is pretty accurate | |
247 | else if((trials > 10) && (trials - 1 > guess) && (guess > 3)) | |
248 | factor = 1.15f; // less accurate but OK. | |
249 | else if(trials < 10) | |
250 | { | |
251 | // pretty inaccurate guess in this area: | |
252 | if(guess > trials / 64) | |
253 | { | |
254 | guess = trials / 4; | |
255 | factor = 2; | |
256 | } | |
257 | else | |
258 | guess = trials / 1024; | |
259 | } | |
260 | else | |
261 | factor = 2; // trials largish, but in far tails. | |
262 | ||
263 | typedef typename Policy::discrete_quantile_type discrete_quantile_type; | |
264 | boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>(); | |
265 | return detail::inverse_discrete_quantile( | |
266 | dist, | |
267 | comp ? q : p, | |
268 | comp, | |
269 | guess, | |
270 | factor, | |
271 | RealType(1), | |
272 | discrete_quantile_type(), | |
273 | max_iter); | |
274 | } // quantile | |
275 | ||
276 | } | |
277 | ||
278 | template <class RealType = double, class Policy = policies::policy<> > | |
279 | class binomial_distribution | |
280 | { | |
281 | public: | |
282 | typedef RealType value_type; | |
283 | typedef Policy policy_type; | |
284 | ||
285 | binomial_distribution(RealType n = 1, RealType p = 0.5) : m_n(n), m_p(p) | |
286 | { // Default n = 1 is the Bernoulli distribution | |
287 | // with equal probability of 'heads' or 'tails. | |
288 | RealType r; | |
289 | binomial_detail::check_dist( | |
290 | "boost::math::binomial_distribution<%1%>::binomial_distribution", | |
291 | m_n, | |
292 | m_p, | |
293 | &r, Policy()); | |
294 | } // binomial_distribution constructor. | |
295 | ||
296 | RealType success_fraction() const | |
297 | { // Probability. | |
298 | return m_p; | |
299 | } | |
300 | RealType trials() const | |
301 | { // Total number of trials. | |
302 | return m_n; | |
303 | } | |
304 | ||
305 | enum interval_type{ | |
306 | clopper_pearson_exact_interval, | |
307 | jeffreys_prior_interval | |
308 | }; | |
309 | ||
310 | // | |
311 | // Estimation of the success fraction parameter. | |
312 | // The best estimate is actually simply successes/trials, | |
313 | // these functions are used | |
314 | // to obtain confidence intervals for the success fraction. | |
315 | // | |
316 | static RealType find_lower_bound_on_p( | |
317 | RealType trials, | |
318 | RealType successes, | |
319 | RealType probability, | |
320 | interval_type t = clopper_pearson_exact_interval) | |
321 | { | |
322 | static const char* function = "boost::math::binomial_distribution<%1%>::find_lower_bound_on_p"; | |
323 | // Error checks: | |
324 | RealType result = 0; | |
325 | if(false == binomial_detail::check_dist_and_k( | |
326 | function, trials, RealType(0), successes, &result, Policy()) | |
327 | && | |
328 | binomial_detail::check_dist_and_prob( | |
329 | function, trials, RealType(0), probability, &result, Policy())) | |
330 | { return result; } | |
331 | ||
332 | if(successes == 0) | |
333 | return 0; | |
334 | ||
335 | // NOTE!!! The Clopper Pearson formula uses "successes" not | |
336 | // "successes+1" as usual to get the lower bound, | |
337 | // see http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm | |
338 | return (t == clopper_pearson_exact_interval) ? ibeta_inv(successes, trials - successes + 1, probability, static_cast<RealType*>(0), Policy()) | |
339 | : ibeta_inv(successes + 0.5f, trials - successes + 0.5f, probability, static_cast<RealType*>(0), Policy()); | |
340 | } | |
341 | static RealType find_upper_bound_on_p( | |
342 | RealType trials, | |
343 | RealType successes, | |
344 | RealType probability, | |
345 | interval_type t = clopper_pearson_exact_interval) | |
346 | { | |
347 | static const char* function = "boost::math::binomial_distribution<%1%>::find_upper_bound_on_p"; | |
348 | // Error checks: | |
349 | RealType result = 0; | |
350 | if(false == binomial_detail::check_dist_and_k( | |
351 | function, trials, RealType(0), successes, &result, Policy()) | |
352 | && | |
353 | binomial_detail::check_dist_and_prob( | |
354 | function, trials, RealType(0), probability, &result, Policy())) | |
355 | { return result; } | |
356 | ||
357 | if(trials == successes) | |
358 | return 1; | |
359 | ||
360 | return (t == clopper_pearson_exact_interval) ? ibetac_inv(successes + 1, trials - successes, probability, static_cast<RealType*>(0), Policy()) | |
361 | : ibetac_inv(successes + 0.5f, trials - successes + 0.5f, probability, static_cast<RealType*>(0), Policy()); | |
362 | } | |
363 | // Estimate number of trials parameter: | |
364 | // | |
365 | // "How many trials do I need to be P% sure of seeing k events?" | |
366 | // or | |
367 | // "How many trials can I have to be P% sure of seeing fewer than k events?" | |
368 | // | |
369 | static RealType find_minimum_number_of_trials( | |
370 | RealType k, // number of events | |
371 | RealType p, // success fraction | |
372 | RealType alpha) // risk level | |
373 | { | |
374 | static const char* function = "boost::math::binomial_distribution<%1%>::find_minimum_number_of_trials"; | |
375 | // Error checks: | |
376 | RealType result = 0; | |
377 | if(false == binomial_detail::check_dist_and_k( | |
378 | function, k, p, k, &result, Policy()) | |
379 | && | |
380 | binomial_detail::check_dist_and_prob( | |
381 | function, k, p, alpha, &result, Policy())) | |
382 | { return result; } | |
383 | ||
384 | result = ibetac_invb(k + 1, p, alpha, Policy()); // returns n - k | |
385 | return result + k; | |
386 | } | |
387 | ||
388 | static RealType find_maximum_number_of_trials( | |
389 | RealType k, // number of events | |
390 | RealType p, // success fraction | |
391 | RealType alpha) // risk level | |
392 | { | |
393 | static const char* function = "boost::math::binomial_distribution<%1%>::find_maximum_number_of_trials"; | |
394 | // Error checks: | |
395 | RealType result = 0; | |
396 | if(false == binomial_detail::check_dist_and_k( | |
397 | function, k, p, k, &result, Policy()) | |
398 | && | |
399 | binomial_detail::check_dist_and_prob( | |
400 | function, k, p, alpha, &result, Policy())) | |
401 | { return result; } | |
402 | ||
403 | result = ibeta_invb(k + 1, p, alpha, Policy()); // returns n - k | |
404 | return result + k; | |
405 | } | |
406 | ||
407 | private: | |
408 | RealType m_n; // Not sure if this shouldn't be an int? | |
409 | RealType m_p; // success_fraction | |
410 | }; // template <class RealType, class Policy> class binomial_distribution | |
411 | ||
412 | typedef binomial_distribution<> binomial; | |
413 | // typedef binomial_distribution<double> binomial; | |
414 | // IS now included since no longer a name clash with function binomial. | |
415 | //typedef binomial_distribution<double> binomial; // Reserved name of type double. | |
416 | ||
417 | template <class RealType, class Policy> | |
418 | const std::pair<RealType, RealType> range(const binomial_distribution<RealType, Policy>& dist) | |
419 | { // Range of permissible values for random variable k. | |
420 | using boost::math::tools::max_value; | |
421 | return std::pair<RealType, RealType>(static_cast<RealType>(0), dist.trials()); | |
422 | } | |
423 | ||
424 | template <class RealType, class Policy> | |
425 | const std::pair<RealType, RealType> support(const binomial_distribution<RealType, Policy>& dist) | |
426 | { // Range of supported values for random variable k. | |
427 | // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. | |
428 | return std::pair<RealType, RealType>(static_cast<RealType>(0), dist.trials()); | |
429 | } | |
430 | ||
431 | template <class RealType, class Policy> | |
432 | inline RealType mean(const binomial_distribution<RealType, Policy>& dist) | |
433 | { // Mean of Binomial distribution = np. | |
434 | return dist.trials() * dist.success_fraction(); | |
435 | } // mean | |
436 | ||
437 | template <class RealType, class Policy> | |
438 | inline RealType variance(const binomial_distribution<RealType, Policy>& dist) | |
439 | { // Variance of Binomial distribution = np(1-p). | |
440 | return dist.trials() * dist.success_fraction() * (1 - dist.success_fraction()); | |
441 | } // variance | |
442 | ||
443 | template <class RealType, class Policy> | |
444 | RealType pdf(const binomial_distribution<RealType, Policy>& dist, const RealType& k) | |
445 | { // Probability Density/Mass Function. | |
446 | BOOST_FPU_EXCEPTION_GUARD | |
447 | ||
448 | BOOST_MATH_STD_USING // for ADL of std functions | |
449 | ||
450 | RealType n = dist.trials(); | |
451 | ||
452 | // Error check: | |
453 | RealType result = 0; // initialization silences some compiler warnings | |
454 | if(false == binomial_detail::check_dist_and_k( | |
455 | "boost::math::pdf(binomial_distribution<%1%> const&, %1%)", | |
456 | n, | |
457 | dist.success_fraction(), | |
458 | k, | |
459 | &result, Policy())) | |
460 | { | |
461 | return result; | |
462 | } | |
463 | ||
464 | // Special cases of success_fraction, regardless of k successes and regardless of n trials. | |
465 | if (dist.success_fraction() == 0) | |
466 | { // probability of zero successes is 1: | |
467 | return static_cast<RealType>(k == 0 ? 1 : 0); | |
468 | } | |
469 | if (dist.success_fraction() == 1) | |
470 | { // probability of n successes is 1: | |
471 | return static_cast<RealType>(k == n ? 1 : 0); | |
472 | } | |
473 | // k argument may be integral, signed, or unsigned, or floating point. | |
474 | // If necessary, it has already been promoted from an integral type. | |
475 | if (n == 0) | |
476 | { | |
477 | return 1; // Probability = 1 = certainty. | |
478 | } | |
479 | if (k == 0) | |
480 | { // binomial coeffic (n 0) = 1, | |
481 | // n ^ 0 = 1 | |
482 | return pow(1 - dist.success_fraction(), n); | |
483 | } | |
484 | if (k == n) | |
485 | { // binomial coeffic (n n) = 1, | |
486 | // n ^ 0 = 1 | |
487 | return pow(dist.success_fraction(), k); // * pow((1 - dist.success_fraction()), (n - k)) = 1 | |
488 | } | |
489 | ||
490 | // Probability of getting exactly k successes | |
491 | // if C(n, k) is the binomial coefficient then: | |
492 | // | |
493 | // f(k; n,p) = C(n, k) * p^k * (1-p)^(n-k) | |
494 | // = (n!/(k!(n-k)!)) * p^k * (1-p)^(n-k) | |
495 | // = (tgamma(n+1) / (tgamma(k+1)*tgamma(n-k+1))) * p^k * (1-p)^(n-k) | |
496 | // = p^k (1-p)^(n-k) / (beta(k+1, n-k+1) * (n+1)) | |
497 | // = ibeta_derivative(k+1, n-k+1, p) / (n+1) | |
498 | // | |
499 | using boost::math::ibeta_derivative; // a, b, x | |
500 | return ibeta_derivative(k+1, n-k+1, dist.success_fraction(), Policy()) / (n+1); | |
501 | ||
502 | ||
503 | ||
504 | template <class RealType, class Policy> | |
505 | inline RealType cdf(const binomial_distribution<RealType, Policy>& dist, const RealType& k) | |
506 | { // Cumulative Distribution Function Binomial. | |
507 | // The random variate k is the number of successes in n trials. | |
508 | // k argument may be integral, signed, or unsigned, or floating point. | |
509 | // If necessary, it has already been promoted from an integral type. | |
510 | ||
511 | // Returns the sum of the terms 0 through k of the Binomial Probability Density/Mass: | |
512 | // | |
513 | // i=k | |
514 | // -- ( n ) i n-i | |
515 | // > | | p (1-p) | |
516 | // -- ( i ) | |
517 | // i=0 | |
518 | ||
519 | // The terms are not summed directly instead | |
520 | // the incomplete beta integral is employed, | |
521 | // according to the formula: | |
522 | // P = I[1-p]( n-k, k+1). | |
523 | // = 1 - I[p](k + 1, n - k) | |
524 | ||
525 | BOOST_MATH_STD_USING // for ADL of std functions | |
526 | ||
527 | RealType n = dist.trials(); | |
528 | RealType p = dist.success_fraction(); | |
529 | ||
530 | // Error check: | |
531 | RealType result = 0; | |
532 | if(false == binomial_detail::check_dist_and_k( | |
533 | "boost::math::cdf(binomial_distribution<%1%> const&, %1%)", | |
534 | n, | |
535 | p, | |
536 | k, | |
537 | &result, Policy())) | |
538 | { | |
539 | return result; | |
540 | } | |
541 | if (k == n) | |
542 | { | |
543 | return 1; | |
544 | } | |
545 | ||
546 | // Special cases, regardless of k. | |
547 | if (p == 0) | |
548 | { // This need explanation: | |
549 | // the pdf is zero for all cases except when k == 0. | |
550 | // For zero p the probability of zero successes is one. | |
551 | // Therefore the cdf is always 1: | |
552 | // the probability of k or *fewer* successes is always 1 | |
553 | // if there are never any successes! | |
554 | return 1; | |
555 | } | |
556 | if (p == 1) | |
557 | { // This is correct but needs explanation: | |
558 | // when k = 1 | |
559 | // all the cdf and pdf values are zero *except* when k == n, | |
560 | // and that case has been handled above already. | |
561 | return 0; | |
562 | } | |
563 | // | |
564 | // P = I[1-p](n - k, k + 1) | |
565 | // = 1 - I[p](k + 1, n - k) | |
566 | // Use of ibetac here prevents cancellation errors in calculating | |
567 | // 1-p if p is very small, perhaps smaller than machine epsilon. | |
568 | // | |
569 | // Note that we do not use a finite sum here, since the incomplete | |
570 | // beta uses a finite sum internally for integer arguments, so | |
571 | // we'll just let it take care of the necessary logic. | |
572 | // | |
573 | return ibetac(k + 1, n - k, p, Policy()); | |
574 | } // binomial cdf | |
575 | ||
576 | template <class RealType, class Policy> | |
577 | inline RealType cdf(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c) | |
578 | { // Complemented Cumulative Distribution Function Binomial. | |
579 | // The random variate k is the number of successes in n trials. | |
580 | // k argument may be integral, signed, or unsigned, or floating point. | |
581 | // If necessary, it has already been promoted from an integral type. | |
582 | ||
583 | // Returns the sum of the terms k+1 through n of the Binomial Probability Density/Mass: | |
584 | // | |
585 | // i=n | |
586 | // -- ( n ) i n-i | |
587 | // > | | p (1-p) | |
588 | // -- ( i ) | |
589 | // i=k+1 | |
590 | ||
591 | // The terms are not summed directly instead | |
592 | // the incomplete beta integral is employed, | |
593 | // according to the formula: | |
594 | // Q = 1 -I[1-p]( n-k, k+1). | |
595 | // = I[p](k + 1, n - k) | |
596 | ||
597 | BOOST_MATH_STD_USING // for ADL of std functions | |
598 | ||
599 | RealType const& k = c.param; | |
600 | binomial_distribution<RealType, Policy> const& dist = c.dist; | |
601 | RealType n = dist.trials(); | |
602 | RealType p = dist.success_fraction(); | |
603 | ||
604 | // Error checks: | |
605 | RealType result = 0; | |
606 | if(false == binomial_detail::check_dist_and_k( | |
607 | "boost::math::cdf(binomial_distribution<%1%> const&, %1%)", | |
608 | n, | |
609 | p, | |
610 | k, | |
611 | &result, Policy())) | |
612 | { | |
613 | return result; | |
614 | } | |
615 | ||
616 | if (k == n) | |
617 | { // Probability of greater than n successes is necessarily zero: | |
618 | return 0; | |
619 | } | |
620 | ||
621 | // Special cases, regardless of k. | |
622 | if (p == 0) | |
623 | { | |
624 | // This need explanation: the pdf is zero for all | |
625 | // cases except when k == 0. For zero p the probability | |
626 | // of zero successes is one. Therefore the cdf is always | |
627 | // 1: the probability of *more than* k successes is always 0 | |
628 | // if there are never any successes! | |
629 | return 0; | |
630 | } | |
631 | if (p == 1) | |
632 | { | |
633 | // This needs explanation, when p = 1 | |
634 | // we always have n successes, so the probability | |
635 | // of more than k successes is 1 as long as k < n. | |
636 | // The k == n case has already been handled above. | |
637 | return 1; | |
638 | } | |
639 | // | |
640 | // Calculate cdf binomial using the incomplete beta function. | |
641 | // Q = 1 -I[1-p](n - k, k + 1) | |
642 | // = I[p](k + 1, n - k) | |
643 | // Use of ibeta here prevents cancellation errors in calculating | |
644 | // 1-p if p is very small, perhaps smaller than machine epsilon. | |
645 | // | |
646 | // Note that we do not use a finite sum here, since the incomplete | |
647 | // beta uses a finite sum internally for integer arguments, so | |
648 | // we'll just let it take care of the necessary logic. | |
649 | // | |
650 | return ibeta(k + 1, n - k, p, Policy()); | |
651 | } // binomial cdf | |
652 | ||
653 | template <class RealType, class Policy> | |
654 | inline RealType quantile(const binomial_distribution<RealType, Policy>& dist, const RealType& p) | |
655 | { | |
656 | return binomial_detail::quantile_imp(dist, p, RealType(1-p), false); | |
657 | } // quantile | |
658 | ||
659 | template <class RealType, class Policy> | |
660 | RealType quantile(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c) | |
661 | { | |
662 | return binomial_detail::quantile_imp(c.dist, RealType(1-c.param), c.param, true); | |
663 | } // quantile | |
664 | ||
665 | template <class RealType, class Policy> | |
666 | inline RealType mode(const binomial_distribution<RealType, Policy>& dist) | |
667 | { | |
668 | BOOST_MATH_STD_USING // ADL of std functions. | |
669 | RealType p = dist.success_fraction(); | |
670 | RealType n = dist.trials(); | |
671 | return floor(p * (n + 1)); | |
672 | } | |
673 | ||
674 | template <class RealType, class Policy> | |
675 | inline RealType median(const binomial_distribution<RealType, Policy>& dist) | |
676 | { // Bounds for the median of the negative binomial distribution | |
677 | // VAN DE VEN R. ; WEBER N. C. ; | |
678 | // Univ. Sydney, school mathematics statistics, Sydney N.S.W. 2006, AUSTRALIE | |
679 | // Metrika (Metrika) ISSN 0026-1335 CODEN MTRKA8 | |
680 | // 1993, vol. 40, no3-4, pp. 185-189 (4 ref.) | |
681 | ||
682 | // Bounds for median and 50 percetage point of binomial and negative binomial distribution | |
683 | // Metrika, ISSN 0026-1335 (Print) 1435-926X (Online) | |
684 | // Volume 41, Number 1 / December, 1994, DOI 10.1007/BF01895303 | |
685 | BOOST_MATH_STD_USING // ADL of std functions. | |
686 | RealType p = dist.success_fraction(); | |
687 | RealType n = dist.trials(); | |
688 | // Wikipedia says one of floor(np) -1, floor (np), floor(np) +1 | |
689 | return floor(p * n); // Chose the middle value. | |
690 | } | |
691 | ||
692 | template <class RealType, class Policy> | |
693 | inline RealType skewness(const binomial_distribution<RealType, Policy>& dist) | |
694 | { | |
695 | BOOST_MATH_STD_USING // ADL of std functions. | |
696 | RealType p = dist.success_fraction(); | |
697 | RealType n = dist.trials(); | |
698 | return (1 - 2 * p) / sqrt(n * p * (1 - p)); | |
699 | } | |
700 | ||
701 | template <class RealType, class Policy> | |
702 | inline RealType kurtosis(const binomial_distribution<RealType, Policy>& dist) | |
703 | { | |
704 | RealType p = dist.success_fraction(); | |
705 | RealType n = dist.trials(); | |
706 | return 3 - 6 / n + 1 / (n * p * (1 - p)); | |
707 | } | |
708 | ||
709 | template <class RealType, class Policy> | |
710 | inline RealType kurtosis_excess(const binomial_distribution<RealType, Policy>& dist) | |
711 | { | |
712 | RealType p = dist.success_fraction(); | |
713 | RealType q = 1 - p; | |
714 | RealType n = dist.trials(); | |
715 | return (1 - 6 * p * q) / (n * p * q); | |
716 | } | |
717 | ||
718 | } // namespace math | |
719 | } // namespace boost | |
720 | ||
721 | // This include must be at the end, *after* the accessors | |
722 | // for this distribution have been defined, in order to | |
723 | // keep compilers that support two-phase lookup happy. | |
724 | #include <boost/math/distributions/detail/derived_accessors.hpp> | |
725 | ||
726 | #endif // BOOST_MATH_SPECIAL_BINOMIAL_HPP | |
727 | ||
728 |