mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
178 lines
7.8 KiB
C++
178 lines
7.8 KiB
C++
// negative_binomial_example2.cpp
|
|
|
|
// Copyright Paul A. Bristow 2007.
|
|
|
|
// 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 examples demonstrating use of the Negative Binomial Distribution.
|
|
// (See other examples for practical applications).
|
|
|
|
#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 distribution: 8 successes (r), 0.25 success fraction = 25% or 1 in 4 successes.
|
|
// negative_binomial_distribution<double> my8dist(8, 0.25);
|
|
negative_binomial my8dist(8, 0.25); // Shorter method using typedef.
|
|
|
|
cout.precision(17); // max_digits10 precision
|
|
cout << "mean(my8dist) = " << mean(my8dist) << endl; // 24
|
|
cout << "my8dist.successes() = " << my8dist.successes() << endl;
|
|
// r th successful trial, after k failures, is r + k th trial.
|
|
cout << "my8dist.success_fraction() = " << my8dist.success_fraction() << endl;
|
|
// success_fraction = failures/successes or k/r = 0.25 or 25%.
|
|
cout << "my8dist.percent success = " << my8dist.success_fraction() * 100 << "%" << endl;
|
|
cout << "cdf(my8dist, 2.) = " << cdf(my8dist, 2.) << endl; // 0.000415802001953125
|
|
cout << "cdf(my8dist, 8.) = " << cdf(my8dist, 8.) << endl; // 0.027129956288263202
|
|
cout << "cdf(complement(my8dist, 8.)) = " << cdf(complement(my8dist, 8.)) << endl; // 0.9728700437117368
|
|
// Check that cdf plus its complement is unity.
|
|
cout << "cdf + complement = " << cdf(my8dist, 8.) + cdf(complement(my8dist, 8.)) << endl; // 1
|
|
// Note: No complement for pdf!
|
|
|
|
// Compare cdf with sum of pdfs.
|
|
double sum = 0.; // of pdfs
|
|
int k = 20;
|
|
for(signed i = 0; i <= k; ++i)
|
|
{
|
|
sum += pdf(my8dist, double(i));
|
|
}
|
|
// Compare with
|
|
double cdf8 = cdf(my8dist, static_cast<double>(k));
|
|
double diff = sum - cdf8;
|
|
|
|
cout << setprecision(17) << "Sum pdfs = " << sum << ' ' // 0.40025683281803698
|
|
<< ", cdf = " << cdf(my8dist, static_cast<double>(k)) // cdf = 0.40025683281803687
|
|
<< ", difference = "
|
|
<< diff/ std::numeric_limits<double>::epsilon() // difference = 0.50000000000000000
|
|
<< " 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 list of values that can be used to plot
|
|
// using Excel, or some other superior graphical display tool.
|
|
|
|
cout << setprecision(17) << showpoint << endl; // max_digits10 precision, including trailing zeros.
|
|
int maxk = static_cast<int>(2. * my8dist.successes() / my8dist.success_fraction());
|
|
// This maxk shows most of the range of interest, probability about 0.0001 to 0.9999.
|
|
cout << " 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(my8dist, static_cast<double>(k)) << ", "
|
|
<< left << setw(25) << cdf(my8dist, static_cast<double>(k))
|
|
<< endl;
|
|
}
|
|
cout << endl;
|
|
|
|
return 0;
|
|
} // int main()
|
|
|
|
/*
|
|
|
|
Output is:
|
|
|
|
negative_binomial distribution - simple example 2
|
|
mean(my8dist) = 24
|
|
my8dist.successes() = 8
|
|
my8dist.success_fraction() = 0.25
|
|
my8dist.percent success = 25%
|
|
cdf(my8dist, 2.) = 0.000415802001953125
|
|
cdf(my8dist, 8.) = 0.027129956288263202
|
|
cdf(complement(my8dist, 8.)) = 0.9728700437117368
|
|
cdf + complement = 1
|
|
Sum pdfs = 0.40025683281803692 , cdf = 0.40025683281803687, difference = 0.25 in epsilon units.
|
|
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
|
|
|
|
|
|
*/
|