mirror of
https://github.com/boostorg/math.git
synced 2026-02-22 03:22:28 +00:00
1002d7d090
Proof reading changes to docs, mostly rewritten neg binomial worked example. Removed dead test. [SVN r3473]
283 lines
10 KiB
Plaintext
283 lines
10 KiB
Plaintext
[section:binom_eg Binomial Distribution Examples]
|
|
|
|
See also the reference documentation for the __binomial_distrib.)
|
|
|
|
[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<>::estimate_lower_bound_on_p` and
|
|
`binomial_distribution<>::estimate_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 occurance" << "= " << 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 Upper\n"
|
|
" Value (%) Limit Limit\n"
|
|
"___________________________________________\n";
|
|
|
|
|
|
And now for the important part - the intervals themselves - for each
|
|
value of /alpha/, we call `estimate_lower_bound_on_p` and
|
|
`estimate_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.
|
|
|
|
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 bounds of interval:
|
|
double l = binomial_distribution<>::estimate_lower_bound_on_p(trials, successes, alpha[i]/2);
|
|
double u = binomial_distribution<>::estimate_upper_bound_on_p(trials, successes, alpha[i]/2);
|
|
// Print Limits:
|
|
cout << fixed << setprecision(5) << setw(15) << right << l;
|
|
cout << fixed << setprecision(5) << setw(15) << right << u << 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 Upper
|
|
Value (%) Limit Limit
|
|
___________________________________________
|
|
50.000 0.12840 0.29588
|
|
75.000 0.09775 0.34633
|
|
90.000 0.07135 0.40103
|
|
95.000 0.05733 0.43661
|
|
99.000 0.03576 0.50661
|
|
99.900 0.01905 0.58632
|
|
99.990 0.01042 0.64997
|
|
99.999 0.00577 0.70216
|
|
''']
|
|
|
|
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]).
|
|
|
|
|
|
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 Upper
|
|
Value (%) Limit Limit
|
|
___________________________________________
|
|
50.000 0.19382 0.20638
|
|
75.000 0.18965 0.21072
|
|
90.000 0.18537 0.21528
|
|
95.000 0.18267 0.21821
|
|
99.000 0.17745 0.22400
|
|
99.900 0.17150 0.23079
|
|
99.990 0.16658 0.23657
|
|
99.999 0.16233 0.24169
|
|
''']
|
|
|
|
Now even when the confidence level is very high, the limits are really
|
|
quite close to the experimentally calculated value of 0.2.
|
|
|
|
[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<>::estimate_number_of_trials`
|
|
can be used to estimate the maximum number of "uses" of that component for some
|
|
probability P.
|
|
|
|
The example program
|
|
[@../../example/binomial_sample_sizes.cpp binomial_sample_sizes.cpp]
|
|
demonstrates its usage. It centres around a routine that prints out
|
|
a table of maximum sample sizes for various probability thresholds:
|
|
|
|
void estimate_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::estimate_number_of_trials(complement(successes, p, alpha[i]));
|
|
t = floor(t);
|
|
// Print Trials:
|
|
cout << fixed << setprecision(5) << setw(15) << right << t << endl;
|
|
}
|
|
|
|
Note that since the value of alpha is the probability of observing
|
|
[*more than /successes/ events], we have to wrap the arguments to
|
|
`estimate_number_of_trials` in a call to complement: remember the binomial
|
|
distribution deals in [*/successes/ or fewer events]. Finally, 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/
|
|
then 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
|
|
(99.9% chance of no failure).
|
|
|
|
[endsect]
|
|
[endsect][/section:binom_eg Binomial Distribution]
|