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