ceph/src/boost/libs/math/example/negative_binomial_example2.cpp

   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 // preceed 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
  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 diference 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 */