// negative_binomial_example2.cpp
// 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)
// Simple example demonstrating use of the Negative Binomial Distribution.
#include <boost/math/distributions/negative_binomial.hpp>
using boost::math::negative_binomial_distribution;
using boost::math::negative_binomial; // typedef
// In a sequence of trials or events
// (Bernoulli, independent, yes or no, succeed or fail)
// with success_fraction probability p,
// negative_binomial is the probability that k or fewer failures
// preceed the r th trial's success.
#include <iostream>
using std::cout;
using std::endl;
using std::setprecision;
using std::showpoint;
using std::setw;
using std::left;
using std::right;
#include <limits>
using std::numeric_limits;
int main()
{
cout << "Negative_binomial distribution - simple example 2" << endl;
// Construct a negative binomial distribution with:
// 8 successes (r), success fraction (p) 0.25 = 25% or 1 in 4 successes.
negative_binomial mynbdist(8, 0.25); // Shorter method using typedef.
// Display (to check) properties of the distribution just constructed.
cout << "mean(mynbdist) = " << mean(mynbdist) << endl; // 24
cout << "mynbdist.successes() = " << mynbdist.successes() << endl; // 8
// r th successful trial, after k failures, is r + k th trial.
cout << "mynbdist.success_fraction() = " << mynbdist.success_fraction() << endl;
// success_fraction = failures/successes or k/r = 0.25 or 25%.
cout << "mynbdist.percent success = " << mynbdist.success_fraction() * 100 << "%" << endl;
// Show as % too.
// Show some cumulative distribution function values for failures k = 2 and 8
cout << "cdf(mynbdist, 2.) = " << cdf(mynbdist, 2.) << endl; // 0.000415802001953125
cout << "cdf(mynbdist, 8.) = " << cdf(mynbdist, 8.) << endl; // 0.027129956288263202
cout << "cdf(complement(mynbdist, 8.)) = " << cdf(complement(mynbdist, 8.)) << endl; // 0.9728700437117368
// Check that cdf plus its complement is unity.
cout << "cdf + complement = " << cdf(mynbdist, 8.) + cdf(complement(mynbdist, 8.)) << endl; // 1
// Note: No complement for pdf!
// Compare cdf with sum of pdfs.
double sum = 0.; // Calculate the sum of all the pdfs,
int k = 20; // for 20 failures
for(signed i = 0; i <= k; ++i)
{
sum += pdf(mynbdist, double(i));
}
// Compare with the cdf
double cdf8 = cdf(mynbdist, static_cast<double>(k));
double diff = sum - cdf8; // Expect the diference to be very small.
cout << setprecision(17) << "Sum pdfs = " << sum << ' ' // sum = 0.40025683281803698
<< ", cdf = " << cdf(mynbdist, static_cast<double>(k)) // cdf = 0.40025683281803687
<< ", difference = " // difference = 0.50000000000000000
<< setprecision(1) << diff/ (std::numeric_limits<double>::epsilon() * sum)
<< " in epsilon units." << endl;
// Note: Use boost::math::tools::epsilon rather than std::numeric_limits
// to cover RealTypes that do not specialize numeric_limits.
//[neg_binomial_example2
// Print a table of values that can be used to plot
// using Excel, or some other superior graphical display tool.
cout.precision(17); // Use max_digits10 precision, the maximum available for a reference table.
cout << showpoint << endl; // include trailing zeros.
// This is a maximum possible precision for the type (here double) to suit a reference table.
int maxk = static_cast<int>(2. * mynbdist.successes() / mynbdist.success_fraction());
// This maxk shows most of the range of interest, probability about 0.0001 to 0.999.
cout << "\n"" k pdf cdf""\n" << endl;
for (int k = 0; k < maxk; k++)
{
cout << right << setprecision(17) << showpoint
<< right << setw(3) << k << ", "
<< left << setw(25) << pdf(mynbdist, static_cast<double>(k))
<< left << setw(25) << cdf(mynbdist, static_cast<double>(k))
<< endl;
}
cout << endl;
//] [/ neg_binomial_example2]
return 0;
} // int main()
/*
Output is:
negative_binomial distribution - simple example 2
mean(mynbdist) = 24
mynbdist.successes() = 8
mynbdist.success_fraction() = 0.25
mynbdist.percent success = 25%
cdf(mynbdist, 2.) = 0.000415802001953125
cdf(mynbdist, 8.) = 0.027129956288263202
cdf(complement(mynbdist, 8.)) = 0.9728700437117368
cdf + complement = 1
Sum pdfs = 0.40025683281803692 , cdf = 0.40025683281803687, difference = 0.25 in epsilon units.
//[neg_binomial_example2_1
k pdf cdf
0, 1.5258789062500000e-005 1.5258789062500003e-005
1, 9.1552734375000000e-005 0.00010681152343750000
2, 0.00030899047851562522 0.00041580200195312500
3, 0.00077247619628906272 0.0011882781982421875
4, 0.0015932321548461918 0.0027815103530883789
5, 0.0028678178787231476 0.0056493282318115234
6, 0.0046602040529251142 0.010309532284736633
7, 0.0069903060793876605 0.017299838364124298
8, 0.0098301179241389001 0.027129956288263202
9, 0.013106823898851871 0.040236780187115073
10, 0.016711200471036140 0.056947980658151209
11, 0.020509200578089786 0.077457181236241013
12, 0.024354675686481652 0.10181185692272265
13, 0.028101548869017230 0.12991340579173993
14, 0.031614242477644432 0.16152764826938440
15, 0.034775666725408917 0.19630331499479325
16, 0.037492515688331451 0.23379583068312471
17, 0.039697957787645101 0.27349378847076977
18, 0.041352039362130305 0.31484582783290005
19, 0.042440250924291580 0.35728607875719176
20, 0.042970754060845245 0.40025683281803687
21, 0.042970754060845225 0.44322758687888220
22, 0.042482450037426581 0.48571003691630876
23, 0.041558918514873783 0.52726895543118257
24, 0.040260202311284021 0.56752915774246648
25, 0.038649794218832620 0.60617895196129912
26, 0.036791631035234917 0.64297058299653398
27, 0.034747651533277427 0.67771823452981139
28, 0.032575923312447595 0.71029415784225891
29, 0.030329307911589130 0.74062346575384819
30, 0.028054609818219924 0.76867807557206813
31, 0.025792141284492545 0.79447021685656061
32, 0.023575629142856460 0.81804584599941710
33, 0.021432390129869489 0.83947823612928651
34, 0.019383705779220189 0.85886194190850684
35, 0.017445335201298231 0.87630727710980494
36, 0.015628112784496322 0.89193538989430121
37, 0.013938587078064250 0.90587397697236549
38, 0.012379666154859701 0.91825364312722524
39, 0.010951243136991251 0.92920488626421649
40, 0.0096507830144735539 0.93885566927869002
41, 0.0084738582566109364 0.94732952753530097
42, 0.0074146259745345548 0.95474415350983555
43, 0.0064662435824429246 0.96121039709227851
44, 0.0056212231142827853 0.96683162020656122
45, 0.0048717266990450708 0.97170334690560634
46, 0.0042098073105878630 0.97591315421619418
47, 0.0036275999165703964 0.97954075413276465
48, 0.0031174686783026818 0.98265822281106729
49, 0.0026721160099737302 0.98533033882104104
50, 0.0022846591885275322 0.98761499800956853
51, 0.0019486798960970148 0.98956367790566557
52, 0.0016582516423517923 0.99122192954801736
53, 0.0014079495076571762 0.99262987905567457
54, 0.0011928461106539983 0.99382272516632852
55, 0.0010084971662802015 0.99483122233260868
56, 0.00085091948404891532 0.99568214181665760
57, 0.00071656377604119542 0.99639870559269883
58, 0.00060228420831048650 0.99700098980100937
59, 0.00050530624256557675 0.99750629604357488
60, 0.00042319397814867202 0.99792949002172360
61, 0.00035381791615708398 0.99828330793788067
62, 0.00029532382517950324 0.99857863176306016
63, 0.00024610318764958566 0.99882473495070978
//] [neg_binomial_example2_1 end of Quickbook]
*/
