{

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 difference 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.

// 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;

return 0;