| [section:binom_eg Binomial Distribution Examples] |
| |
| See also the reference documentation for the __binomial_distrib. |
| |
| [section:binomial_coinflip_example Binomial Coin-Flipping Example] |
| |
| [import ../../../example/binomial_coinflip_example.cpp] |
| [binomial_coinflip_example1] |
| |
| See [@../../../example/binomial_coinflip_example.cpp binomial_coinflip_example.cpp] |
| for full source code, the program output looks like this: |
| |
| [binomial_coinflip_example_output] |
| |
| [endsect] [/section:binomial_coinflip_example Binomial coinflip example] |
| |
| [section:binomial_quiz_example Binomial Quiz Example] |
| |
| [import ../../../example/binomial_quiz_example.cpp] |
| [binomial_quiz_example1] |
| [binomial_quiz_example2] |
| [discrete_quantile_real] |
| |
| See [@../../../example/binomial_quiz_example.cpp binomial_quiz_example.cpp] |
| for full source code and output. |
| |
| [endsect] [/section:binomial_coinflip_quiz Binomial Coin-Flipping example] |
| |
| [section:binom_conf Calculating Confidence Limits on the Frequency of Occurrence for a Binomial Distribution] |
| |
| Imagine you have a process that follows a binomial distribution: for each |
| trial conducted, an event either occurs or does it does not, referred |
| to as "successes" and "failures". If, by experiment, you want to measure the |
| frequency with which successes occur, the best estimate is given simply |
| by /k/ \/ /N/, for /k/ successes out of /N/ trials. However our confidence in that |
| estimate will be shaped by how many trials were conducted, and how many successes |
| were observed. The static member functions |
| `binomial_distribution<>::find_lower_bound_on_p` and |
| `binomial_distribution<>::find_upper_bound_on_p` allow you to calculate |
| the confidence intervals for your estimate of the occurrence frequency. |
| |
| The sample program [@../../../example/binomial_confidence_limits.cpp |
| binomial_confidence_limits.cpp] illustrates their use. It begins by defining |
| a procedure that will print a table of confidence limits for various degrees |
| of certainty: |
| |
| #include <iostream> |
| #include <iomanip> |
| #include <boost/math/distributions/binomial.hpp> |
| |
| void confidence_limits_on_frequency(unsigned trials, unsigned successes) |
| { |
| // |
| // trials = Total number of trials. |
| // successes = Total number of observed successes. |
| // |
| // Calculate confidence limits for an observed |
| // frequency of occurrence that follows a binomial |
| // distribution. |
| // |
| using namespace std; |
| using namespace boost::math; |
| |
| // Print out general info: |
| cout << |
| "___________________________________________\n" |
| "2-Sided Confidence Limits For Success Ratio\n" |
| "___________________________________________\n\n"; |
| cout << setprecision(7); |
| cout << setw(40) << left << "Number of Observations" << "= " << trials << "\n"; |
| cout << setw(40) << left << "Number of successes" << "= " << successes << "\n"; |
| cout << setw(40) << left << "Sample frequency of occurrence" << "= " << double(successes) / trials << "\n"; |
| |
| The procedure now defines a table of significance levels: these are the |
| probabilities that the true occurrence frequency lies outside the calculated |
| interval: |
| |
| double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 }; |
| |
| Some pretty printing of the table header follows: |
| |
| cout << "\n\n" |
| "_______________________________________________________________________\n" |
| "Confidence Lower CP Upper CP Lower JP Upper JP\n" |
| " Value (%) Limit Limit Limit Limit\n" |
| "_______________________________________________________________________\n"; |
| |
| |
| And now for the important part - the intervals themselves - for each |
| value of /alpha/, we call `find_lower_bound_on_p` and |
| `find_lower_upper_on_p` to obtain lower and upper bounds |
| respectively. Note that since we are calculating a two-sided interval, |
| we must divide the value of alpha in two. |
| |
| Please note that calculating two separate /single sided bounds/, each with risk |
| level [alpha][space]is not the same thing as calculating a two sided interval. |
| Had we calculate two single-sided intervals each with a risk |
| that the true value is outside the interval of [alpha], then: |
| |
| * The risk that it is less than the lower bound is [alpha]. |
| |
| and |
| |
| * The risk that it is greater than the upper bound is also [alpha]. |
| |
| So the risk it is outside *upper or lower bound*, is *twice* alpha, and the |
| probability that it is inside the bounds is therefore not nearly as high as |
| one might have thought. This is why [alpha]/2 must be used in |
| the calculations below. |
| |
| In contrast, had we been calculating a |
| single-sided interval, for example: ['"Calculate a lower bound so that we are P% |
| sure that the true occurrence frequency is greater than some value"] |
| then we would *not* have divided by two. |
| |
| Finally note that `binomial_distribution` provides a choice of two |
| methods for the calculation, we print out the results from both |
| methods in this example: |
| |
| for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i) |
| { |
| // Confidence value: |
| cout << fixed << setprecision(3) << setw(10) << right << 100 * (1-alpha[i]); |
| // Calculate Clopper Pearson bounds: |
| double l = binomial_distribution<>::find_lower_bound_on_p( |
| trials, successes, alpha[i]/2); |
| double u = binomial_distribution<>::find_upper_bound_on_p( |
| trials, successes, alpha[i]/2); |
| // Print Clopper Pearson Limits: |
| cout << fixed << setprecision(5) << setw(15) << right << l; |
| cout << fixed << setprecision(5) << setw(15) << right << u; |
| // Calculate Jeffreys Prior Bounds: |
| l = binomial_distribution<>::find_lower_bound_on_p( |
| trials, successes, alpha[i]/2, |
| binomial_distribution<>::jeffreys_prior_interval); |
| u = binomial_distribution<>::find_upper_bound_on_p( |
| trials, successes, alpha[i]/2, |
| binomial_distribution<>::jeffreys_prior_interval); |
| // Print Jeffreys Prior Limits: |
| cout << fixed << setprecision(5) << setw(15) << right << l; |
| cout << fixed << setprecision(5) << setw(15) << right << u << std::endl; |
| } |
| cout << endl; |
| } |
| |
| And that's all there is to it. Let's see some sample output for a 2 in 10 |
| success ratio, first for 20 trials: |
| |
| [pre'''___________________________________________ |
| 2-Sided Confidence Limits For Success Ratio |
| ___________________________________________ |
| |
| Number of Observations = 20 |
| Number of successes = 4 |
| Sample frequency of occurrence = 0.2 |
| |
| |
| _______________________________________________________________________ |
| Confidence Lower CP Upper CP Lower JP Upper JP |
| Value (%) Limit Limit Limit Limit |
| _______________________________________________________________________ |
| 50.000 0.12840 0.29588 0.14974 0.26916 |
| 75.000 0.09775 0.34633 0.11653 0.31861 |
| 90.000 0.07135 0.40103 0.08734 0.37274 |
| 95.000 0.05733 0.43661 0.07152 0.40823 |
| 99.000 0.03576 0.50661 0.04655 0.47859 |
| 99.900 0.01905 0.58632 0.02634 0.55960 |
| 99.990 0.01042 0.64997 0.01530 0.62495 |
| 99.999 0.00577 0.70216 0.00901 0.67897 |
| '''] |
| |
| As you can see, even at the 95% confidence level the bounds are |
| really quite wide (this example is chosen to be easily compared to the one |
| in the __handbook |
| [@http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm |
| here]). Note also that the Clopper-Pearson calculation method (CP above) |
| produces quite noticeably more pessimistic estimates than the Jeffreys Prior |
| method (JP above). |
| |
| |
| Compare that with the program output for |
| 2000 trials: |
| |
| [pre'''___________________________________________ |
| 2-Sided Confidence Limits For Success Ratio |
| ___________________________________________ |
| |
| Number of Observations = 2000 |
| Number of successes = 400 |
| Sample frequency of occurrence = 0.2000000 |
| |
| |
| _______________________________________________________________________ |
| Confidence Lower CP Upper CP Lower JP Upper JP |
| Value (%) Limit Limit Limit Limit |
| _______________________________________________________________________ |
| 50.000 0.19382 0.20638 0.19406 0.20613 |
| 75.000 0.18965 0.21072 0.18990 0.21047 |
| 90.000 0.18537 0.21528 0.18561 0.21503 |
| 95.000 0.18267 0.21821 0.18291 0.21796 |
| 99.000 0.17745 0.22400 0.17769 0.22374 |
| 99.900 0.17150 0.23079 0.17173 0.23053 |
| 99.990 0.16658 0.23657 0.16681 0.23631 |
| 99.999 0.16233 0.24169 0.16256 0.24143 |
| '''] |
| |
| Now even when the confidence level is very high, the limits are really |
| quite close to the experimentally calculated value of 0.2. Furthermore |
| the difference between the two calculation methods is now really quite small. |
| |
| [endsect] |
| |
| [section:binom_size_eg Estimating Sample Sizes for a Binomial Distribution.] |
| |
| Imagine you have a critical component that you know will fail in 1 in |
| N "uses" (for some suitable definition of "use"). You may want to schedule |
| routine replacement of the component so that its chance of failure between |
| routine replacements is less than P%. If the failures follow a binomial |
| distribution (each time the component is "used" it either fails or does not) |
| then the static member function `binomial_distibution<>::find_maximum_number_of_trials` |
| can be used to estimate the maximum number of "uses" of that component for some |
| acceptable risk level /alpha/. |
| |
| The example program |
| [@../../../example/binomial_sample_sizes.cpp binomial_sample_sizes.cpp] |
| demonstrates its usage. It centres on a routine that prints out |
| a table of maximum sample sizes for various probability thresholds: |
| |
| void find_max_sample_size( |
| double p, // success ratio. |
| unsigned successes) // Total number of observed successes permitted. |
| { |
| |
| The routine then declares a table of probability thresholds: these are the |
| maximum acceptable probability that /successes/ or fewer events will be |
| observed. In our example, /successes/ will be always zero, since we want |
| no component failures, but in other situations non-zero values may well |
| make sense. |
| |
| double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 }; |
| |
| Much of the rest of the program is pretty-printing, the important part |
| is in the calculation of maximum number of permitted trials for each |
| value of alpha: |
| |
| for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i) |
| { |
| // Confidence value: |
| cout << fixed << setprecision(3) << setw(10) << right << 100 * (1-alpha[i]); |
| // calculate trials: |
| double t = binomial::find_maximum_number_of_trials( |
| successes, p, alpha[i]); |
| t = floor(t); |
| // Print Trials: |
| cout << fixed << setprecision(5) << setw(15) << right << t << endl; |
| } |
| |
| Note that since we're |
| calculating the maximum number of trials permitted, we'll err on the safe |
| side and take the floor of the result. Had we been calculating the |
| /minimum/ number of trials required to observe a certain number of /successes/ |
| using `find_minimum_number_of_trials` we would have taken the ceiling instead. |
| |
| We'll finish off by looking at some sample output, firstly for |
| a 1 in 1000 chance of component failure with each use: |
| |
| [pre |
| '''________________________ |
| Maximum Number of Trials |
| ________________________ |
| |
| Success ratio = 0.001 |
| Maximum Number of "successes" permitted = 0 |
| |
| |
| ____________________________ |
| Confidence Max Number |
| Value (%) Of Trials |
| ____________________________ |
| 50.000 692 |
| 75.000 287 |
| 90.000 105 |
| 95.000 51 |
| 99.000 10 |
| 99.900 0 |
| 99.990 0 |
| 99.999 0''' |
| ] |
| |
| So 51 "uses" of the component would yield a 95% chance that no |
| component failures would be observed. |
| |
| Compare that with a 1 in 1 million chance of component failure: |
| |
| [pre''' |
| ________________________ |
| Maximum Number of Trials |
| ________________________ |
| |
| Success ratio = 0.0000010 |
| Maximum Number of "successes" permitted = 0 |
| |
| |
| ____________________________ |
| Confidence Max Number |
| Value (%) Of Trials |
| ____________________________ |
| 50.000 693146 |
| 75.000 287681 |
| 90.000 105360 |
| 95.000 51293 |
| 99.000 10050 |
| 99.900 1000 |
| 99.990 100 |
| 99.999 10''' |
| ] |
| |
| In this case, even 1000 uses of the component would still yield a |
| less than 1 in 1000 chance of observing a component failure |
| (i.e. a 99.9% chance of no failure). |
| |
| [endsect] [/section:binom_size_eg Estimating Sample Sizes for a Binomial Distribution.] |
| |
| [endsect][/section:binom_eg Binomial Distribution] |
| |
| [/ |
| Copyright 2006 John Maddock and Paul A. Bristow. |
| Distributed under the Boost Software License, Version 1.0. |
| (See accompanying file LICENSE_1_0.txt or copy at |
| http://www.boost.org/LICENSE_1_0.txt). |
| ] |
| |