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 // precede 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
main()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 difference 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 */
183