mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
Naive example giving listing suitable for plotting etc.
[SVN r3302]
This commit is contained in:
Paul A. Bristow
204
example/negative_binomial_example2.cpp
Normal file
204
example/negative_binomial_example2.cpp
Normal file
@@ -0,0 +1,204 @@
|
||||
// negative_binomial_example2.cpp
|
||||
|
||||
// Copyright Paul A. Bristow 2006.
|
||||
|
||||
// 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 for Negative Binomial Distribution.
|
||||
|
||||
#define BOOST_MATH_THROW_ON_DOMAIN_ERROR
|
||||
|
||||
#ifdef _MSC_VER
|
||||
# pragma warning(disable: 4127) // conditional expression is constant.
|
||||
# pragma warning(disable: 4100) // unreferenced formal parameter.
|
||||
# pragma warning(disable: 4512) // assignment operator could not be generated.
|
||||
#endif
|
||||
|
||||
#include <boost/math/distributions/negative_binomial.hpp> // for negative_binomial_distribution
|
||||
using boost::math::negative_binomial_distribution;
|
||||
|
||||
// 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 <boost/math/special_functions/gamma.hpp>
|
||||
using boost::math::lgamma; // log gamma
|
||||
|
||||
//#include <boost/math/concepts/real_concept.hpp> // for real_concept
|
||||
//using ::boost::math::concepts::real_concept;
|
||||
|
||||
//#include <boost/test/included/test_exec_monitor.hpp> // for test_main
|
||||
//#include <boost/test/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE
|
||||
|
||||
#include <iostream>
|
||||
using std::cout;
|
||||
using std::endl;
|
||||
using std::setprecision;
|
||||
using std::showpoint;
|
||||
using std::setw;
|
||||
using std::left;
|
||||
#include <limits>
|
||||
using std::numeric_limits;
|
||||
|
||||
int main(int, char* [])
|
||||
{
|
||||
#ifdef BOOST_MATH_THROW_ON_DOMAIN_ERROR
|
||||
cout << "BOOST_MATH_THROW_ON_DOMAIN_ERROR" << " is defined to throw on domain error." << endl;
|
||||
#else
|
||||
cout << "BOOST_MATH_THROW_ON_DOMAIN_ERROR" << " is NOT defined, so NO throw on domain error." << endl;
|
||||
#endif
|
||||
cout << "negative_binomial distribution - simples examples" << endl;
|
||||
using namespace boost::math;
|
||||
cout << setprecision(17) << showpoint << endl; // max_digits10 precision including trailing zeros.
|
||||
|
||||
negative_binomial_distribution<double> my8dist(8., 0.25);
|
||||
// 8 successes (r), 0.25 success fraction = 25% or 1 in 4 successes.
|
||||
// Note: double values (matching the distribution definition) avoid the need for any casting.
|
||||
|
||||
cout << "mean(my8dist) = " << mean(my8dist) << endl; //
|
||||
cout << "my8dist.successes() = " << my8dist.successes() << endl; // r th trial is successful.
|
||||
cout << "my8dist.success_fraction() = " << my8dist.success_fraction() << endl; //
|
||||
cout << "cdf(my8dist, 2.) = " << cdf(my8dist, 2.) << endl; // 4.1580200195313E-4
|
||||
|
||||
cout << "cdf(my8dist, 8.) = " << cdf(my8dist, 8.) << endl;
|
||||
cout << "cdf(complement(my8dist, 8.)) = " << cdf(complement(my8dist, 8.)) << endl;
|
||||
cout << "cdf + complement = " << cdf(my8dist, 8.) + cdf(complement(my8dist, 8.)) << endl;
|
||||
// Expect unity.
|
||||
// cout << "cdf(complement(my8dist, 8.)) = " << pdf(complement(my8dist, 8.)) << endl;
|
||||
// Note: No complement for pdf!
|
||||
double sum = 0.;
|
||||
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) << showpoint << "Sum pdfs = " << sum <<' ' // 0.40025683281803714
|
||||
<< ", cdf = " << cdf(my8dist, static_cast<double>(k))
|
||||
<< ", difference = "
|
||||
<< diff/ std::numeric_limits<double>::epsilon()
|
||||
<< " in epsilon units." << endl;
|
||||
// 0.40025683281803681
|
||||
//double tol = boost::math::tools::epsilon<double>() * 100 * 10; // 100 for %, so 10 eps
|
||||
// Use boost::math::tools::epsilon rather than std::numeric_limits to cover
|
||||
// RealTypes that do not specialize numeric_limits.
|
||||
//BOOST_CHECK_CLOSE(sum, cdf(my8dist, static_cast<double>(20)), tol);
|
||||
|
||||
|
||||
// Print a list of values that can be used to plot
|
||||
// using Excel, or some superior graphical display tool.
|
||||
|
||||
int maxk = static_cast<int>(2. * my8dist.successes() / my8dist.success_fraction());
|
||||
// This shows most of the range of interest.
|
||||
for (int k = 0; k < maxk; k++)
|
||||
{
|
||||
cout << setprecision(17) << showpoint
|
||||
<< left << setw(3) << k << ", "
|
||||
<< left << setw(23) << pdf(my8dist, static_cast<double>(k)) << ", "
|
||||
<< left << setw(23) << cdf(my8dist, static_cast<double>(k))
|
||||
<< endl;
|
||||
}
|
||||
|
||||
|
||||
return 0;
|
||||
} // int main(int, char* [])
|
||||
|
||||
/*
|
||||
|
||||
Output is
|
||||
|
||||
------ Build started: Project: negative_binomial_example2, Configuration: Debug Win32 ------
|
||||
Compiling...
|
||||
negative_binomial_example2.cpp
|
||||
Linking...
|
||||
Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\negative_binomial_example2.exe"
|
||||
BOOST_MATH_THROW_ON_DOMAIN_ERROR is defined to throw on domain error.
|
||||
negative_binomial distribution - simples examples
|
||||
mean(my8dist) = 24.000000000000000
|
||||
my8dist.successes() = 8.0000000000000000
|
||||
my8dist.success_fraction() = 0.25000000000000000
|
||||
cdf(my8dist, 2.) = 0.00041580200195312495
|
||||
cdf(my8dist, 8.) = 0.027129956288263198
|
||||
cdf(complement(my8dist, 8.)) = 0.97287004371173680
|
||||
cdf + complement = 1.0000000000000000
|
||||
Sum pdfs = 0.40025683281803714 , cdf = 0.40025683281803681, difference = 1.5000000000000000 in epsilon units.
|
||||
0 , 1.5258789062499998e-005, 1.5258789062500003e-005
|
||||
1 , 9.1552734374999959e-005, 0.00010681152343750000
|
||||
2 , 0.00030899047851562522 , 0.00041580200195312495
|
||||
3 , 0.00077247619628906239 , 0.0011882781982421873
|
||||
4 , 0.0015932321548461931 , 0.0027815103530883785
|
||||
5 , 0.0028678178787231463 , 0.0056493282318115226
|
||||
6 , 0.0046602040529251142 , 0.010309532284736632
|
||||
7 , 0.0069903060793876605 , 0.017299838364124295
|
||||
8 , 0.0098301179241389071 , 0.027129956288263198
|
||||
9 , 0.013106823898851858 , 0.040236780187115066
|
||||
10 , 0.016711200471036133 , 0.056947980658151202
|
||||
11 , 0.020509200578089803 , 0.077457181236240999
|
||||
12 , 0.024354675686481628 , 0.10181185692272264
|
||||
13 , 0.028101548869017345 , 0.12991340579173991
|
||||
14 , 0.031614242477644432 , 0.16152764826938434
|
||||
15 , 0.034775666725408945 , 0.19630331499479323
|
||||
16 , 0.037492515688331500 , 0.23379583068312468
|
||||
17 , 0.039697957787645122 , 0.27349378847076972
|
||||
18 , 0.041352039362130291 , 0.31484582783289999
|
||||
19 , 0.042440250924291600 , 0.35728607875719159
|
||||
20 , 0.042970754060845266 , 0.40025683281803681
|
||||
21 , 0.042970754060845252 , 0.44322758687888208
|
||||
22 , 0.042482450037426567 , 0.48571003691630860
|
||||
23 , 0.041558918514873790 , 0.52726895543118235
|
||||
24 , 0.040260202311284035 , 0.56752915774246648
|
||||
25 , 0.038649794218832613 , 0.60617895196129901
|
||||
26 , 0.036791631035234952 , 0.64297058299653465
|
||||
27 , 0.034747651533277447 , 0.67771823452981150
|
||||
28 , 0.032575923312447616 , 0.71029415784225836
|
||||
29 , 0.030329307911589130 , 0.74062346575384863
|
||||
30 , 0.028054609818219937 , 0.76867807557206813
|
||||
31 , 0.025792141284492535 , 0.79447021685656094
|
||||
32 , 0.023575629142856481 , 0.81804584599941688
|
||||
33 , 0.021432390129869510 , 0.83947823612928674
|
||||
34 , 0.019383705779220179 , 0.85886194190850684
|
||||
35 , 0.017445335201298227 , 0.87630727710980472
|
||||
36 , 0.015628112784496315 , 0.89193538989430132
|
||||
37 , 0.013938587078064250 , 0.90587397697236538
|
||||
38 , 0.012379666154859706 , 0.91825364312722502
|
||||
39 , 0.010951243136991251 , 0.92920488626421649
|
||||
40 , 0.0096507830144735348 , 0.93885566927868991
|
||||
41 , 0.0084738582566109191 , 0.94732952753530109
|
||||
42 , 0.0074146259745345999 , 0.95474415350983566
|
||||
43 , 0.0064662435824429246 , 0.96121039709227840
|
||||
44 , 0.0056212231142827853 , 0.96683162020656122
|
||||
45 , 0.0048717266990450708 , 0.97170334690560634
|
||||
46 , 0.0042098073105878457 , 0.97591315421619418
|
||||
47 , 0.0036275999165703964 , 0.97954075413276454
|
||||
48 , 0.0031174686783026818 , 0.98265822281106729
|
||||
49 , 0.0026721160099737293 , 0.98533033882104104
|
||||
50 , 0.0022846591885275322 , 0.98761499800956853
|
||||
51 , 0.0019486798960970076 , 0.98956367790566557
|
||||
52 , 0.0016582516423517852 , 0.99122192954801736
|
||||
53 , 0.0014079495076571762 , 0.99262987905567457
|
||||
54 , 0.0011928461106539983 , 0.99382272516632852
|
||||
55 , 0.0010084971662801955 , 0.99483122233260879
|
||||
56 , 0.00085091948404891532 , 0.99568214181665760
|
||||
57 , 0.00071656377604119553 , 0.99639870559269883
|
||||
58 , 0.00060228420831048911 , 0.99700098980100937
|
||||
59 , 0.00050530624256557902 , 0.99750629604357488
|
||||
60 , 0.00042319397814867321 , 0.99792949002172360
|
||||
61 , 0.00035381791615708398 , 0.99828330793788067
|
||||
62 , 0.00029532382517950264 , 0.99857863176306016
|
||||
63 , 0.00024610318764958615 , 0.99882473495070978
|
||||
Build Time 0:03
|
||||
Build log was saved at "file://i:\boost-06-05-03-1300\libs\math\test\Math_test\neative_binomial_example2\Debug\BuildLog.htm"
|
||||
negative_binomial_example2 - 0 error(s), 0 warning(s)
|
||||
========== Build: 1 succeeded, 0 failed, 0 up-to-date, 0 skipped ==========
|
||||
|
||||
|
||||
*/
|
||||
Reference in New Issue
Block a user