1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
|
// neg_binomial_confidence_limits.cpp
// Copyright John Maddock 2006
// Copyright Paul A. Bristow 2007, 2010
// Use, modification and distribution are subject to 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)
// Caution: this file contains quickbook markup as well as code
// and comments, don't change any of the special comment markups!
//[neg_binomial_confidence_limits
/*`
First we need some includes to access the negative binomial distribution
(and some basic std output of course).
*/
#include <boost/math/distributions/negative_binomial.hpp>
using boost::math::negative_binomial;
#include <iostream>
using std::cout; using std::endl;
#include <iomanip>
using std::setprecision;
using std::setw; using std::left; using std::fixed; using std::right;
/*`
First define 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 };
/*`
Confidence value as % is (1 - alpha) * 100, so alpha 0.05 == 95% confidence
that the true occurrence frequency lies *inside* the calculated interval.
We need a function to calculate and print confidence limits
for an observed frequency of occurrence
that follows a negative binomial distribution.
*/
void confidence_limits_on_frequency(unsigned trials, unsigned successes)
{
// trials = Total number of trials.
// successes = Total number of observed successes.
// failures = trials - successes.
// success_fraction = successes /trials.
// Print out general info:
cout <<
"______________________________________________\n"
"2-Sided Confidence Limits For Success Fraction\n"
"______________________________________________\n\n";
cout << setprecision(7);
cout << setw(40) << left << "Number of trials" << " = " << trials << "\n";
cout << setw(40) << left << "Number of successes" << " = " << successes << "\n";
cout << setw(40) << left << "Number of failures" << " = " << trials - successes << "\n";
cout << setw(40) << left << "Observed frequency of occurrence" << " = " << double(successes) / trials << "\n";
// Print table header:
cout << "\n\n"
"___________________________________________\n"
"Confidence Lower Upper\n"
" Value (%) Limit Limit\n"
"___________________________________________\n";
/*`
And now for the important part - the bounds themselves.
For each value of /alpha/, we call `find_lower_bound_on_p` and
`find_upper_bound_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. 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.
*/
// Now print out the upper and lower limits for the alpha table values.
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:
double lower = negative_binomial::find_lower_bound_on_p(trials, successes, alpha[i]/2);
double upper = negative_binomial::find_upper_bound_on_p(trials, successes, alpha[i]/2);
// Print limits:
cout << fixed << setprecision(5) << setw(15) << right << lower;
cout << fixed << setprecision(5) << setw(15) << right << upper << endl;
}
cout << endl;
} // void confidence_limits_on_frequency(unsigned trials, unsigned successes)
/*`
And then call confidence_limits_on_frequency with increasing numbers of trials,
but always the same success fraction 0.1, or 1 in 10.
*/
int main()
{
confidence_limits_on_frequency(20, 2); // 20 trials, 2 successes, 2 in 20, = 1 in 10 = 0.1 success fraction.
confidence_limits_on_frequency(200, 20); // More trials, but same 0.1 success fraction.
confidence_limits_on_frequency(2000, 200); // Many more trials, but same 0.1 success fraction.
return 0;
} // int main()
//] [/negative_binomial_confidence_limits_eg end of Quickbook in C++ markup]
/*
______________________________________________
2-Sided Confidence Limits For Success Fraction
______________________________________________
Number of trials = 20
Number of successes = 2
Number of failures = 18
Observed frequency of occurrence = 0.1
___________________________________________
Confidence Lower Upper
Value (%) Limit Limit
___________________________________________
50.000 0.04812 0.13554
75.000 0.03078 0.17727
90.000 0.01807 0.22637
95.000 0.01235 0.26028
99.000 0.00530 0.33111
99.900 0.00164 0.41802
99.990 0.00051 0.49202
99.999 0.00016 0.55574
______________________________________________
2-Sided Confidence Limits For Success Fraction
______________________________________________
Number of trials = 200
Number of successes = 20
Number of failures = 180
Observed frequency of occurrence = 0.1000000
___________________________________________
Confidence Lower Upper
Value (%) Limit Limit
___________________________________________
50.000 0.08462 0.11350
75.000 0.07580 0.12469
90.000 0.06726 0.13695
95.000 0.06216 0.14508
99.000 0.05293 0.16170
99.900 0.04343 0.18212
99.990 0.03641 0.20017
99.999 0.03095 0.21664
______________________________________________
2-Sided Confidence Limits For Success Fraction
______________________________________________
Number of trials = 2000
Number of successes = 200
Number of failures = 1800
Observed frequency of occurrence = 0.1000000
___________________________________________
Confidence Lower Upper
Value (%) Limit Limit
___________________________________________
50.000 0.09536 0.10445
75.000 0.09228 0.10776
90.000 0.08916 0.11125
95.000 0.08720 0.11352
99.000 0.08344 0.11802
99.900 0.07921 0.12336
99.990 0.07577 0.12795
99.999 0.07282 0.13206
*/
|
