]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // negative_binomial_example2.cpp |
2 | ||
3 | // Copyright Paul A. Bristow 2007, 2010. | |
4 | ||
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 | // Simple example demonstrating use of the Negative Binomial Distribution. | |
11 | ||
12 | #include <boost/math/distributions/negative_binomial.hpp> | |
13 | using boost::math::negative_binomial_distribution; | |
14 | using boost::math::negative_binomial; // typedef | |
15 | ||
16 | // In a sequence of trials or events | |
17 | // (Bernoulli, independent, yes or no, succeed or fail) | |
18 | // with success_fraction probability p, | |
19 | // negative_binomial is the probability that k or fewer failures | |
20 | // preceed the r th trial's success. | |
21 | ||
22 | #include <iostream> | |
23 | using std::cout; | |
24 | using std::endl; | |
25 | using std::setprecision; | |
26 | using std::showpoint; | |
27 | using std::setw; | |
28 | using std::left; | |
29 | using std::right; | |
30 | #include <limits> | |
31 | using std::numeric_limits; | |
32 | ||
33 | int main() | |
34 | { | |
35 | cout << "Negative_binomial distribution - simple example 2" << endl; | |
36 | // Construct a negative binomial distribution with: | |
37 | // 8 successes (r), success fraction (p) 0.25 = 25% or 1 in 4 successes. | |
38 | negative_binomial mynbdist(8, 0.25); // Shorter method using typedef. | |
39 | ||
40 | // Display (to check) properties of the distribution just constructed. | |
41 | cout << "mean(mynbdist) = " << mean(mynbdist) << endl; // 24 | |
42 | cout << "mynbdist.successes() = " << mynbdist.successes() << endl; // 8 | |
43 | // r th successful trial, after k failures, is r + k th trial. | |
44 | cout << "mynbdist.success_fraction() = " << mynbdist.success_fraction() << endl; | |
45 | // success_fraction = failures/successes or k/r = 0.25 or 25%. | |
46 | cout << "mynbdist.percent success = " << mynbdist.success_fraction() * 100 << "%" << endl; | |
47 | // Show as % too. | |
48 | // Show some cumulative distribution function values for failures k = 2 and 8 | |
49 | cout << "cdf(mynbdist, 2.) = " << cdf(mynbdist, 2.) << endl; // 0.000415802001953125 | |
50 | cout << "cdf(mynbdist, 8.) = " << cdf(mynbdist, 8.) << endl; // 0.027129956288263202 | |
51 | cout << "cdf(complement(mynbdist, 8.)) = " << cdf(complement(mynbdist, 8.)) << endl; // 0.9728700437117368 | |
52 | // Check that cdf plus its complement is unity. | |
53 | cout << "cdf + complement = " << cdf(mynbdist, 8.) + cdf(complement(mynbdist, 8.)) << endl; // 1 | |
54 | // Note: No complement for pdf! | |
55 | ||
56 | // Compare cdf with sum of pdfs. | |
57 | double sum = 0.; // Calculate the sum of all the pdfs, | |
58 | int k = 20; // for 20 failures | |
59 | for(signed i = 0; i <= k; ++i) | |
60 | { | |
61 | sum += pdf(mynbdist, double(i)); | |
62 | } | |
63 | // Compare with the cdf | |
64 | double cdf8 = cdf(mynbdist, static_cast<double>(k)); | |
65 | double diff = sum - cdf8; // Expect the diference to be very small. | |
66 | cout << setprecision(17) << "Sum pdfs = " << sum << ' ' // sum = 0.40025683281803698 | |
67 | << ", cdf = " << cdf(mynbdist, static_cast<double>(k)) // cdf = 0.40025683281803687 | |
68 | << ", difference = " // difference = 0.50000000000000000 | |
69 | << setprecision(1) << diff/ (std::numeric_limits<double>::epsilon() * sum) | |
70 | << " in epsilon units." << endl; | |
71 | ||
72 | // Note: Use boost::math::tools::epsilon rather than std::numeric_limits | |
73 | // to cover RealTypes that do not specialize numeric_limits. | |
74 | ||
75 | //[neg_binomial_example2 | |
76 | ||
77 | // Print a table of values that can be used to plot | |
78 | // using Excel, or some other superior graphical display tool. | |
79 | ||
80 | cout.precision(17); // Use max_digits10 precision, the maximum available for a reference table. | |
81 | cout << showpoint << endl; // include trailing zeros. | |
82 | // This is a maximum possible precision for the type (here double) to suit a reference table. | |
83 | int maxk = static_cast<int>(2. * mynbdist.successes() / mynbdist.success_fraction()); | |
84 | // This maxk shows most of the range of interest, probability about 0.0001 to 0.999. | |
85 | cout << "\n"" k pdf cdf""\n" << endl; | |
86 | for (int k = 0; k < maxk; k++) | |
87 | { | |
88 | cout << right << setprecision(17) << showpoint | |
89 | << right << setw(3) << k << ", " | |
90 | << left << setw(25) << pdf(mynbdist, static_cast<double>(k)) | |
91 | << left << setw(25) << cdf(mynbdist, static_cast<double>(k)) | |
92 | << endl; | |
93 | } | |
94 | cout << endl; | |
95 | //] [/ neg_binomial_example2] | |
96 | return 0; | |
97 | } // int main() | |
98 | ||
99 | /* | |
100 | ||
101 | Output is: | |
102 | ||
103 | negative_binomial distribution - simple example 2 | |
104 | mean(mynbdist) = 24 | |
105 | mynbdist.successes() = 8 | |
106 | mynbdist.success_fraction() = 0.25 | |
107 | mynbdist.percent success = 25% | |
108 | cdf(mynbdist, 2.) = 0.000415802001953125 | |
109 | cdf(mynbdist, 8.) = 0.027129956288263202 | |
110 | cdf(complement(mynbdist, 8.)) = 0.9728700437117368 | |
111 | cdf + complement = 1 | |
112 | Sum pdfs = 0.40025683281803692 , cdf = 0.40025683281803687, difference = 0.25 in epsilon units. | |
113 | ||
114 | //[neg_binomial_example2_1 | |
115 | k pdf cdf | |
116 | 0, 1.5258789062500000e-005 1.5258789062500003e-005 | |
117 | 1, 9.1552734375000000e-005 0.00010681152343750000 | |
118 | 2, 0.00030899047851562522 0.00041580200195312500 | |
119 | 3, 0.00077247619628906272 0.0011882781982421875 | |
120 | 4, 0.0015932321548461918 0.0027815103530883789 | |
121 | 5, 0.0028678178787231476 0.0056493282318115234 | |
122 | 6, 0.0046602040529251142 0.010309532284736633 | |
123 | 7, 0.0069903060793876605 0.017299838364124298 | |
124 | 8, 0.0098301179241389001 0.027129956288263202 | |
125 | 9, 0.013106823898851871 0.040236780187115073 | |
126 | 10, 0.016711200471036140 0.056947980658151209 | |
127 | 11, 0.020509200578089786 0.077457181236241013 | |
128 | 12, 0.024354675686481652 0.10181185692272265 | |
129 | 13, 0.028101548869017230 0.12991340579173993 | |
130 | 14, 0.031614242477644432 0.16152764826938440 | |
131 | 15, 0.034775666725408917 0.19630331499479325 | |
132 | 16, 0.037492515688331451 0.23379583068312471 | |
133 | 17, 0.039697957787645101 0.27349378847076977 | |
134 | 18, 0.041352039362130305 0.31484582783290005 | |
135 | 19, 0.042440250924291580 0.35728607875719176 | |
136 | 20, 0.042970754060845245 0.40025683281803687 | |
137 | 21, 0.042970754060845225 0.44322758687888220 | |
138 | 22, 0.042482450037426581 0.48571003691630876 | |
139 | 23, 0.041558918514873783 0.52726895543118257 | |
140 | 24, 0.040260202311284021 0.56752915774246648 | |
141 | 25, 0.038649794218832620 0.60617895196129912 | |
142 | 26, 0.036791631035234917 0.64297058299653398 | |
143 | 27, 0.034747651533277427 0.67771823452981139 | |
144 | 28, 0.032575923312447595 0.71029415784225891 | |
145 | 29, 0.030329307911589130 0.74062346575384819 | |
146 | 30, 0.028054609818219924 0.76867807557206813 | |
147 | 31, 0.025792141284492545 0.79447021685656061 | |
148 | 32, 0.023575629142856460 0.81804584599941710 | |
149 | 33, 0.021432390129869489 0.83947823612928651 | |
150 | 34, 0.019383705779220189 0.85886194190850684 | |
151 | 35, 0.017445335201298231 0.87630727710980494 | |
152 | 36, 0.015628112784496322 0.89193538989430121 | |
153 | 37, 0.013938587078064250 0.90587397697236549 | |
154 | 38, 0.012379666154859701 0.91825364312722524 | |
155 | 39, 0.010951243136991251 0.92920488626421649 | |
156 | 40, 0.0096507830144735539 0.93885566927869002 | |
157 | 41, 0.0084738582566109364 0.94732952753530097 | |
158 | 42, 0.0074146259745345548 0.95474415350983555 | |
159 | 43, 0.0064662435824429246 0.96121039709227851 | |
160 | 44, 0.0056212231142827853 0.96683162020656122 | |
161 | 45, 0.0048717266990450708 0.97170334690560634 | |
162 | 46, 0.0042098073105878630 0.97591315421619418 | |
163 | 47, 0.0036275999165703964 0.97954075413276465 | |
164 | 48, 0.0031174686783026818 0.98265822281106729 | |
165 | 49, 0.0026721160099737302 0.98533033882104104 | |
166 | 50, 0.0022846591885275322 0.98761499800956853 | |
167 | 51, 0.0019486798960970148 0.98956367790566557 | |
168 | 52, 0.0016582516423517923 0.99122192954801736 | |
169 | 53, 0.0014079495076571762 0.99262987905567457 | |
170 | 54, 0.0011928461106539983 0.99382272516632852 | |
171 | 55, 0.0010084971662802015 0.99483122233260868 | |
172 | 56, 0.00085091948404891532 0.99568214181665760 | |
173 | 57, 0.00071656377604119542 0.99639870559269883 | |
174 | 58, 0.00060228420831048650 0.99700098980100937 | |
175 | 59, 0.00050530624256557675 0.99750629604357488 | |
176 | 60, 0.00042319397814867202 0.99792949002172360 | |
177 | 61, 0.00035381791615708398 0.99828330793788067 | |
178 | 62, 0.00029532382517950324 0.99857863176306016 | |
179 | 63, 0.00024610318764958566 0.99882473495070978 | |
180 | //] [neg_binomial_example2_1 end of Quickbook] | |
181 | ||
182 | */ |