libs/math/test/test_poisson.cpp

   1 // test_poisson.cpp
   2
   3 // Copyright Paul A. Bristow 2007.
   4 // Copyright John Maddock 2006.
   5
   6 // Use, modification and distribution are subject to the
   7 // Boost Software License, Version 1.0.
   8 // (See accompanying file LICENSE_1_0.txt
   9 // or copy at http://www.boost.org/LICENSE_1_0.txt)
  10
  11 // Basic sanity test for Poisson Cumulative Distribution Function.
  12
  13 #define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
  14
  15 #if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
  16 #  define TEST_FLOAT
  17 #  define TEST_DOUBLE
  18 #  define TEST_LDOUBLE
  19 #  define TEST_REAL_CONCEPT
  20 #endif
  21
  22 #ifdef _MSC_VER
  23 #  pragma warning(disable: 4127) // conditional expression is constant.
  24 #endif
  25
  26 #define BOOST_TEST_MAIN
  27 #include <boost/test/unit_test.hpp> // Boost.Test
  28 #include <boost/test/tools/floating_point_comparison.hpp>
  29
  30 #include <boost/math/concepts/real_concept.hpp> // for real_concept
  31 #include <boost/math/distributions/poisson.hpp>
  32     using boost::math::poisson_distribution;
  33 #include <boost/math/tools/test.hpp> // for real_concept
  34
  35 #include <boost/math/special_functions/gamma.hpp> // for (incomplete) gamma.
  36 //   using boost::math::qamma_Q;
  37 #include "table_type.hpp"
  38 #include "test_out_of_range.hpp"
  39
  40 #include <iostream>
  41    using std::cout;
  42    using std::endl;
  43    using std::setprecision;
  44    using std::showpoint;
  45    using std::ios;
  46 #include <limits>
  47   using std::numeric_limits;
  48
  49 template <class RealType> // Any floating-point type RealType.
  50 void test_spots(RealType)
  51 {
  52   // Basic sanity checks, tolerance is about numeric_limits<RealType>::digits10 decimal places,
  53    // guaranteed for type RealType, eg 6 for float, 15 for double,
  54    // expressed as a percentage (so -2) for BOOST_CHECK_CLOSE,
  55
  56    int decdigits = numeric_limits<RealType>::digits10;
  57   // May eb >15 for 80 and 128-bit FP typtes.
  58   if (decdigits <= 0)
  59   { // decdigits is not defined, for example real concept,
  60     // so assume precision of most test data is double (for example, MathCAD).
  61      decdigits = numeric_limits<double>::digits10; // == 15 for 64-bit
  62   }
  63   if (decdigits > 15 ) // numeric_limits<double>::digits10)
  64   { // 15 is the accuracy of the MathCAD test data.
  65     decdigits = 15; // numeric_limits<double>::digits10;
  66   }
  67
  68    decdigits -= 1; // Perhaps allow some decimal digit(s) margin of numerical error.
  69    RealType tolerance = static_cast<RealType>(std::pow(10., static_cast<double>(2-decdigits))); // 1e-6 (-2 so as %)
  70    tolerance *= 2; // Allow some bit(s) small margin (2 means + or - 1 bit) of numerical error.
  71    // Typically 2e-13% = 2e-15 as fraction for double.
  72
  73    // Sources of spot test values:
  74
  75   // Many be some combinations for which the result is 'exact',
  76   // or at least is good to 40 decimal digits.
  77    // 40 decimal digits includes 128-bit significand User Defined Floating-Point types,
  78
  79    // Best source of accurate values is:
  80    // Mathworld online calculator (40 decimal digits precision, suitable for up to 128-bit significands)
  81    // http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=GammaRegularized
  82    // GammaRegularized is same as gamma incomplete, gamma or gamma_q(a, x) or Q(a, z).
  83
  84   // http://documents.wolfram.com/calculationcenter/v2/Functions/ListsMatrices/Statistics/PoissonDistribution.html
  85
  86   // MathCAD defines ppois(k, lambda== mean) as k integer, k >=0.
  87   // ppois(0, 5) =  6.73794699908547e-3
  88   // ppois(1, 5) = 0.040427681994513;
  89   // ppois(10, 10) = 5.830397501929850E-001
  90   // ppois(10, 1) = 9.999999899522340E-001
  91   // ppois(5,5) = 0.615960654833065
  92
  93   // qpois returns inverse Poission distribution, that is the smallest (floor) k so that ppois(k, lambda) >= p
  94   // p is real number, real mean lambda > 0
  95   // k is approximately the integer for which probability(X <= k) = p
  96   // when random variable X has the Poisson distribution with parameters lambda.
  97   // Uses discrete bisection.
  98   // qpois(6.73794699908547e-3, 5) = 1
  99   // qpois(0.040427681994513, 5) =
 100
 101   // Test Poisson with spot values from MathCAD 'known good'.
 102
 103   using boost::math::poisson_distribution;
 104   using  ::boost::math::poisson;
 105   using  ::boost::math::cdf;
 106   using  ::boost::math::pdf;
 107
 108    // Check that bad arguments throw.
 109    BOOST_MATH_CHECK_THROW(
 110    cdf(poisson_distribution<RealType>(static_cast<RealType>(0)), // mean zero is bad.
 111       static_cast<RealType>(0)),  // even for a good k.
 112       std::domain_error); // Expected error to be thrown.
 113
 114     BOOST_MATH_CHECK_THROW(
 115    cdf(poisson_distribution<RealType>(static_cast<RealType>(-1)), // mean negative is bad.
 116       static_cast<RealType>(0)),
 117       std::domain_error);
 118
 119    BOOST_MATH_CHECK_THROW(
 120    cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unit OK,
 121       static_cast<RealType>(-1)),  // but negative events is bad.
 122       std::domain_error);
 123
 124   BOOST_MATH_CHECK_THROW(
 125      cdf(poisson_distribution<RealType>(static_cast<RealType>(0)), // mean zero is bad.
 126       static_cast<RealType>(99999)),  // for any k events.
 127       std::domain_error);
 128
 129   BOOST_MATH_CHECK_THROW(
 130      cdf(poisson_distribution<RealType>(static_cast<RealType>(0)), // mean zero is bad.
 131       static_cast<RealType>(99999)),  // for any k events.
 132       std::domain_error);
 133
 134   BOOST_MATH_CHECK_THROW(
 135      quantile(poisson_distribution<RealType>(static_cast<RealType>(0)), // mean zero.
 136       static_cast<RealType>(0.5)),  // probability OK.
 137       std::domain_error);
 138
 139   BOOST_MATH_CHECK_THROW(
 140      quantile(poisson_distribution<RealType>(static_cast<RealType>(-1)),
 141       static_cast<RealType>(-1)),  // bad probability.
 142       std::domain_error);
 143
 144   BOOST_MATH_CHECK_THROW(
 145      quantile(poisson_distribution<RealType>(static_cast<RealType>(1)),
 146       static_cast<RealType>(-1)),  // bad probability.
 147       std::domain_error);
 148
 149   BOOST_MATH_CHECK_THROW(
 150      quantile(poisson_distribution<RealType>(static_cast<RealType>(1)),
 151       static_cast<RealType>(1)),  // bad probability.
 152       std::overflow_error);
 153
 154   BOOST_MATH_CHECK_THROW(
 155      quantile(complement(poisson_distribution<RealType>(static_cast<RealType>(1)),
 156       static_cast<RealType>(0))),  // bad probability.
 157       std::overflow_error);
 158
 159   BOOST_CHECK_EQUAL(
 160      quantile(poisson_distribution<RealType>(static_cast<RealType>(1)),
 161       static_cast<RealType>(0)),  // bad probability.
 162       0);
 163
 164   BOOST_CHECK_EQUAL(
 165      quantile(complement(poisson_distribution<RealType>(static_cast<RealType>(1)),
 166       static_cast<RealType>(1))),  // bad probability.
 167       0);
 168
 169   // Check some test values.
 170
 171   BOOST_CHECK_CLOSE( // mode
 172      mode(poisson_distribution<RealType>(static_cast<RealType>(4))), // mode = mean = 4.
 173       static_cast<RealType>(4), // mode.
 174          tolerance);
 175
 176   //BOOST_CHECK_CLOSE( // mode
 177   //   median(poisson_distribution<RealType>(static_cast<RealType>(4))), // mode = mean = 4.
 178   //    static_cast<RealType>(4), // mode.
 179       //   tolerance);
 180   poisson_distribution<RealType> dist4(static_cast<RealType>(40));
 181
 182   BOOST_CHECK_CLOSE( // median
 183      median(dist4), // mode = mean = 4. median = 40.328333333333333
 184       quantile(dist4, static_cast<RealType>(0.5)), // 39.332839138842637
 185          tolerance);
 186
 187   // PDF
 188   BOOST_CHECK_CLOSE(
 189      pdf(poisson_distribution<RealType>(static_cast<RealType>(4)), // mean 4.
 190       static_cast<RealType>(0)),
 191       static_cast<RealType>(1.831563888873410E-002), // probability.
 192          tolerance);
 193
 194   BOOST_CHECK_CLOSE(
 195      pdf(poisson_distribution<RealType>(static_cast<RealType>(4)), // mean 4.
 196       static_cast<RealType>(2)),
 197       static_cast<RealType>(1.465251111098740E-001), // probability.
 198          tolerance);
 199
 200   BOOST_CHECK_CLOSE(
 201      pdf(poisson_distribution<RealType>(static_cast<RealType>(20)), // mean big.
 202       static_cast<RealType>(1)),   //  k small
 203       static_cast<RealType>(4.122307244877130E-008), // probability.
 204          tolerance);
 205
 206   BOOST_CHECK_CLOSE(
 207      pdf(poisson_distribution<RealType>(static_cast<RealType>(4)), // mean 4.
 208       static_cast<RealType>(20)),   //  K>> mean
 209       static_cast<RealType>(8.277463646553730E-009), // probability.
 210          tolerance);
 211
 212   // CDF
 213   BOOST_CHECK_CLOSE(
 214      cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
 215       static_cast<RealType>(0)),  // zero k events.
 216       static_cast<RealType>(3.678794411714420E-1), // probability.
 217          tolerance);
 218
 219   BOOST_CHECK_CLOSE(
 220      cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
 221       static_cast<RealType>(1)),  // one k event.
 222       static_cast<RealType>(7.357588823428830E-1), // probability.
 223          tolerance);
 224
 225   BOOST_CHECK_CLOSE(
 226      cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
 227       static_cast<RealType>(2)),  // two k events.
 228       static_cast<RealType>(9.196986029286060E-1), // probability.
 229          tolerance);
 230
 231   BOOST_CHECK_CLOSE(
 232      cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
 233       static_cast<RealType>(10)),  // two k events.
 234       static_cast<RealType>(9.999999899522340E-1), // probability.
 235          tolerance);
 236
 237   BOOST_CHECK_CLOSE(
 238      cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
 239       static_cast<RealType>(15)),  // two k events.
 240       static_cast<RealType>(9.999999999999810E-1), // probability.
 241          tolerance);
 242
 243   BOOST_CHECK_CLOSE(
 244      cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
 245       static_cast<RealType>(16)),  // two k events.
 246       static_cast<RealType>(9.999999999999990E-1), // probability.
 247          tolerance);
 248
 249   BOOST_CHECK_CLOSE(
 250      cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
 251       static_cast<RealType>(17)),  // two k events.
 252       static_cast<RealType>(1.), // probability unity for double.
 253          tolerance);
 254
 255   BOOST_CHECK_CLOSE(
 256      cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
 257       static_cast<RealType>(33)),  // k events at limit for float unchecked_factorial table.
 258       static_cast<RealType>(1.), // probability.
 259          tolerance);
 260
 261   BOOST_CHECK_CLOSE(
 262      cdf(poisson_distribution<RealType>(static_cast<RealType>(100)), // mean 100.
 263       static_cast<RealType>(33)),  // k events at limit for float unchecked_factorial table.
 264       static_cast<RealType>(6.328271240363390E-15), // probability is tiny.
 265          tolerance * static_cast<RealType>(2e11)); // 6.3495253382825722e-015 MathCAD
 266       // Note that there two tiny probability are much more different.
 267
 268    BOOST_CHECK_CLOSE(
 269      cdf(poisson_distribution<RealType>(static_cast<RealType>(100)), // mean 100.
 270       static_cast<RealType>(34)),  // k events at limit for float unchecked_factorial table.
 271       static_cast<RealType>(1.898481372109020E-14), // probability is tiny.
 272          tolerance*static_cast<RealType>(2e11)); //         1.8984813721090199e-014 MathCAD
 273
 274
 275  BOOST_CHECK_CLOSE(
 276      cdf(poisson_distribution<RealType>(static_cast<RealType>(33)), // mean = k
 277       static_cast<RealType>(33)),  // k events above limit for float unchecked_factorial table.
 278       static_cast<RealType>(5.461191812386560E-1), // probability.
 279          tolerance);
 280
 281  BOOST_CHECK_CLOSE(
 282      cdf(poisson_distribution<RealType>(static_cast<RealType>(33)), // mean = k-1
 283       static_cast<RealType>(34)),  // k events above limit for float unchecked_factorial table.
 284       static_cast<RealType>(6.133535681502950E-1), // probability.
 285          tolerance);
 286
 287  BOOST_CHECK_CLOSE(
 288      cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
 289       static_cast<RealType>(34)),  // k events above limit for float unchecked_factorial table.
 290       static_cast<RealType>(1.), // probability.
 291          tolerance);
 292
 293   BOOST_CHECK_CLOSE(
 294      cdf(poisson_distribution<RealType>(static_cast<RealType>(5.)), // mean
 295       static_cast<RealType>(5)),  // k events.
 296       static_cast<RealType>(0.615960654833065), // probability.
 297          tolerance);
 298   BOOST_CHECK_CLOSE(
 299      cdf(poisson_distribution<RealType>(static_cast<RealType>(5.)), // mean
 300       static_cast<RealType>(1)),  // k events.
 301       static_cast<RealType>(0.040427681994512805), // probability.
 302          tolerance);
 303
 304   BOOST_CHECK_CLOSE(
 305      cdf(poisson_distribution<RealType>(static_cast<RealType>(5.)), // mean
 306       static_cast<RealType>(0)),  // k events (uses special case formula, not gamma).
 307       static_cast<RealType>(0.006737946999085467), // probability.
 308          tolerance);
 309
 310   BOOST_CHECK_CLOSE(
 311      cdf(poisson_distribution<RealType>(static_cast<RealType>(1.)), // mean
 312       static_cast<RealType>(0)),  // k events (uses special case formula, not gamma).
 313       static_cast<RealType>(0.36787944117144233), // probability.
 314          tolerance);
 315
 316   BOOST_CHECK_CLOSE(
 317      cdf(poisson_distribution<RealType>(static_cast<RealType>(10.)), // mean
 318       static_cast<RealType>(10)),  // k events.
 319       static_cast<RealType>(0.5830397501929856), // probability.
 320          tolerance);
 321
 322   BOOST_CHECK_CLOSE(
 323      cdf(poisson_distribution<RealType>(static_cast<RealType>(4.)), // mean
 324       static_cast<RealType>(5)),  // k events.
 325       static_cast<RealType>(0.785130387030406), // probability.
 326          tolerance);
 327
 328   // complement CDF
 329   BOOST_CHECK_CLOSE( // Complement CDF
 330      cdf(complement(poisson_distribution<RealType>(static_cast<RealType>(4.)), // mean
 331       static_cast<RealType>(5))),  // k events.
 332       static_cast<RealType>(1 - 0.785130387030406), // probability.
 333          tolerance);
 334
 335   BOOST_CHECK_CLOSE( // Complement CDF
 336      cdf(complement(poisson_distribution<RealType>(static_cast<RealType>(4.)), // mean
 337       static_cast<RealType>(0))),  // Zero k events (uses special case formula, not gamma).
 338       static_cast<RealType>(0.98168436111126578), // probability.
 339          tolerance);
 340   BOOST_CHECK_CLOSE( // Complement CDF
 341      cdf(complement(poisson_distribution<RealType>(static_cast<RealType>(1.)), // mean
 342       static_cast<RealType>(0))),  // Zero k events (uses special case formula, not gamma).
 343       static_cast<RealType>(0.63212055882855767), // probability.
 344          tolerance);
 345
 346   // Example where k is bigger than max_factorial (>34 for float)
 347   // (therefore using log gamma so perhaps less accurate).
 348   BOOST_CHECK_CLOSE(
 349      cdf(poisson_distribution<RealType>(static_cast<RealType>(40.)), // mean
 350       static_cast<RealType>(40)),  // k events.
 351       static_cast<RealType>(0.5419181783625430), // probability.
 352          tolerance);
 353
 354    // Quantile & complement.
 355   BOOST_CHECK_CLOSE(
 356     boost::math::quantile(
 357          poisson_distribution<RealType>(5),  // mean.
 358          static_cast<RealType>(0.615960654833065)),  //  probability.
 359          static_cast<RealType>(5.), // Expect k = 5
 360          tolerance/5); //
 361
 362   // EQUAL is too optimistic - fails [5.0000000000000124 != 5]
 363   // BOOST_CHECK_EQUAL(boost::math::quantile( //
 364   //       poisson_distribution<RealType>(5.),  // mean.
 365   //       static_cast<RealType>(0.615960654833065)),  //  probability.
 366   //       static_cast<RealType>(5.)); // Expect k = 5 events.
 367
 368   BOOST_CHECK_CLOSE(boost::math::quantile(
 369          poisson_distribution<RealType>(4),  // mean.
 370          static_cast<RealType>(0.785130387030406)),  //  probability.
 371          static_cast<RealType>(5.), // Expect k = 5 events.
 372          tolerance/5);
 373
 374   // Check on quantile of other examples of inverse of cdf.
 375   BOOST_CHECK_CLOSE(
 376      cdf(poisson_distribution<RealType>(static_cast<RealType>(10.)), // mean
 377       static_cast<RealType>(10)),  // k events.
 378       static_cast<RealType>(0.5830397501929856), // probability.
 379          tolerance);
 380
 381   BOOST_CHECK_CLOSE(boost::math::quantile( // inverse of cdf above.
 382          poisson_distribution<RealType>(10.),  // mean.
 383          static_cast<RealType>(0.5830397501929856)),  //  probability.
 384          static_cast<RealType>(10.), // Expect k = 10 events.
 385          tolerance/5);
 386
 387
 388   BOOST_CHECK_CLOSE(
 389      cdf(poisson_distribution<RealType>(static_cast<RealType>(4.)), // mean
 390       static_cast<RealType>(5)),  // k events.
 391       static_cast<RealType>(0.785130387030406), // probability.
 392          tolerance);
 393
 394   BOOST_CHECK_CLOSE(boost::math::quantile( // inverse of cdf above.
 395          poisson_distribution<RealType>(4.),  // mean.
 396          static_cast<RealType>(0.785130387030406)),  //  probability.
 397          static_cast<RealType>(5.), // Expect k = 10 events.
 398          tolerance/5);
 399
 400
 401
 402   //BOOST_CHECK_CLOSE(boost::math::quantile(
 403   //       poisson_distribution<RealType>(5),  // mean.
 404   //       static_cast<RealType>(0.785130387030406)),  //  probability.
 405   //        // 6.1882832344329559 result but MathCAD givest smallest integer ppois(k, mean) >= prob
 406   //       static_cast<RealType>(6.), // Expect k = 6 events.
 407   //       tolerance/5);
 408
 409   //BOOST_CHECK_CLOSE(boost::math::quantile(
 410   //       poisson_distribution<RealType>(5),  // mean.
 411   //       static_cast<RealType>(0.77)),  //  probability.
 412   //        // 6.1882832344329559 result but MathCAD givest smallest integer ppois(k, mean) >= prob
 413   //       static_cast<RealType>(7.), // Expect k = 6 events.
 414   //       tolerance/5);
 415
 416   //BOOST_CHECK_CLOSE(boost::math::quantile(
 417   //       poisson_distribution<RealType>(5),  // mean.
 418   //       static_cast<RealType>(0.75)),  //  probability.
 419   //        // 6.1882832344329559 result but MathCAD givest smallest integer ppois(k, mean) >= prob
 420   //       static_cast<RealType>(6.), // Expect k = 6 events.
 421   //       tolerance/5);
 422
 423   BOOST_CHECK_CLOSE(
 424     boost::math::quantile(
 425          complement(
 426            poisson_distribution<RealType>(4),
 427            static_cast<RealType>(1 - 0.785130387030406))),  // complement.
 428            static_cast<RealType>(5), // Expect k = 5 events.
 429          tolerance/5);
 430
 431   BOOST_CHECK_EQUAL(boost::math::quantile( // Check case when probability < cdf(0) (== pdf(0))
 432          poisson_distribution<RealType>(1),  // mean is small, so cdf and pdf(0) are about 0.35.
 433          static_cast<RealType>(0.0001)),  //  probability < cdf(0).
 434          static_cast<RealType>(0)); // Expect k = 0 events exactly.
 435
 436   BOOST_CHECK_EQUAL(
 437     boost::math::quantile(
 438          complement(
 439            poisson_distribution<RealType>(1),
 440            static_cast<RealType>(0.9999))),  // complement, so 1-probability < cdf(0)
 441            static_cast<RealType>(0)); // Expect k = 0 events exactly.
 442
 443   //
 444   // Test quantile policies against test data:
 445   //
 446 #define T RealType
 447 #include "poisson_quantile.ipp"
 448
 449   for(unsigned i = 0; i < poisson_quantile_data.size(); ++i)
 450   {
 451      using namespace boost::math::policies;
 452      typedef policy<discrete_quantile<boost::math::policies::real> > P1;
 453      typedef policy<discrete_quantile<integer_round_down> > P2;
 454      typedef policy<discrete_quantile<integer_round_up> > P3;
 455      typedef policy<discrete_quantile<integer_round_outwards> > P4;
 456      typedef policy<discrete_quantile<integer_round_inwards> > P5;
 457      typedef policy<discrete_quantile<integer_round_nearest> > P6;
 458      RealType tol = boost::math::tools::epsilon<RealType>() * 20;
 459      if(!boost::is_floating_point<RealType>::value)
 460         tol *= 7;
 461      //
 462      // Check full real value first:
 463      //
 464      poisson_distribution<RealType, P1> p1(poisson_quantile_data[i][0]);
 465      RealType x = quantile(p1, poisson_quantile_data[i][1]);
 466      BOOST_CHECK_CLOSE_FRACTION(x, poisson_quantile_data[i][2], tol);
 467      x = quantile(complement(p1, poisson_quantile_data[i][1]));
 468      BOOST_CHECK_CLOSE_FRACTION(x, poisson_quantile_data[i][3], tol * 3);
 469      //
 470      // Now with round down to integer:
 471      //
 472      poisson_distribution<RealType, P2> p2(poisson_quantile_data[i][0]);
 473      x = quantile(p2, poisson_quantile_data[i][1]);
 474      BOOST_CHECK_EQUAL(x, floor(poisson_quantile_data[i][2]));
 475      x = quantile(complement(p2, poisson_quantile_data[i][1]));
 476      BOOST_CHECK_EQUAL(x, floor(poisson_quantile_data[i][3]));
 477      //
 478      // Now with round up to integer:
 479      //
 480      poisson_distribution<RealType, P3> p3(poisson_quantile_data[i][0]);
 481      x = quantile(p3, poisson_quantile_data[i][1]);
 482      BOOST_CHECK_EQUAL(x, ceil(poisson_quantile_data[i][2]));
 483      x = quantile(complement(p3, poisson_quantile_data[i][1]));
 484      BOOST_CHECK_EQUAL(x, ceil(poisson_quantile_data[i][3]));
 485      //
 486      // Now with round to integer "outside":
 487      //
 488      poisson_distribution<RealType, P4> p4(poisson_quantile_data[i][0]);
 489      x = quantile(p4, poisson_quantile_data[i][1]);
 490      BOOST_CHECK_EQUAL(x, poisson_quantile_data[i][1] < 0.5f ? floor(poisson_quantile_data[i][2]) : ceil(poisson_quantile_data[i][2]));
 491      x = quantile(complement(p4, poisson_quantile_data[i][1]));
 492      BOOST_CHECK_EQUAL(x, poisson_quantile_data[i][1] < 0.5f ? ceil(poisson_quantile_data[i][3]) : floor(poisson_quantile_data[i][3]));
 493      //
 494      // Now with round to integer "inside":
 495      //
 496      poisson_distribution<RealType, P5> p5(poisson_quantile_data[i][0]);
 497      x = quantile(p5, poisson_quantile_data[i][1]);
 498      BOOST_CHECK_EQUAL(x, poisson_quantile_data[i][1] < 0.5f ? ceil(poisson_quantile_data[i][2]) : floor(poisson_quantile_data[i][2]));
 499      x = quantile(complement(p5, poisson_quantile_data[i][1]));
 500      BOOST_CHECK_EQUAL(x, poisson_quantile_data[i][1] < 0.5f ? floor(poisson_quantile_data[i][3]) : ceil(poisson_quantile_data[i][3]));
 501      //
 502      // Now with round to nearest integer:
 503      //
 504      poisson_distribution<RealType, P6> p6(poisson_quantile_data[i][0]);
 505      x = quantile(p6, poisson_quantile_data[i][1]);
 506      BOOST_CHECK_EQUAL(x, floor(poisson_quantile_data[i][2] + 0.5f));
 507      x = quantile(complement(p6, poisson_quantile_data[i][1]));
 508      BOOST_CHECK_EQUAL(x, floor(poisson_quantile_data[i][3] + 0.5f));
 509   }
 510    check_out_of_range<poisson_distribution<RealType> >(1);
 511 } // template <class RealType>void test_spots(RealType)
 512
 513 //
 514
 515 BOOST_AUTO_TEST_CASE( test_main )
 516 {
 517   // Check that can construct normal distribution using the two convenience methods:
 518   using namespace boost::math;
 519   poisson myp1(2); // Using typedef
 520    poisson_distribution<> myp2(2); // Using default RealType double.
 521
 522    // Basic sanity-check spot values.
 523
 524   // Some plain double examples & tests:
 525   cout.precision(17); // double max_digits10
 526   cout.setf(ios::showpoint);
 527
 528   poisson mypoisson(4.); // // mean = 4, default FP type is double.
 529   cout << "mean(mypoisson, 4.) == " << mean(mypoisson) << endl;
 530   cout << "mean(mypoisson, 0.) == " << mean(mypoisson) << endl;
 531   cout << "cdf(mypoisson, 2.) == " << cdf(mypoisson, 2.) << endl;
 532   cout << "pdf(mypoisson, 2.) == " << pdf(mypoisson, 2.) << endl;
 533
 534   // poisson mydudpoisson(0.);
 535   // throws (if BOOST_MATH_DOMAIN_ERROR_POLICY == throw_on_error).
 536
 537 #ifndef BOOST_NO_EXCEPTIONS
 538   BOOST_MATH_CHECK_THROW(poisson mydudpoisson(-1), std::domain_error);// Mean must be > 0.
 539   BOOST_MATH_CHECK_THROW(poisson mydudpoisson(-1), std::logic_error);// Mean must be > 0.
 540 #else
 541   BOOST_MATH_CHECK_THROW(poisson(-1), std::domain_error);// Mean must be > 0.
 542   BOOST_MATH_CHECK_THROW(poisson(-1), std::logic_error);// Mean must be > 0.
 543 #endif
 544   // Passes the check because logic_error is a parent????
 545   // BOOST_MATH_CHECK_THROW(poisson mydudpoisson(-1), std::overflow_error); // fails the check
 546   // because overflow_error is unrelated - except from std::exception
 547   BOOST_MATH_CHECK_THROW(cdf(mypoisson, -1), std::domain_error); // k must be >= 0
 548
 549   BOOST_CHECK_EQUAL(mean(mypoisson), 4.);
 550   BOOST_CHECK_CLOSE(
 551   pdf(mypoisson, 2.),  // k events = 2.
 552     1.465251111098740E-001, // probability.
 553       5e-13);
 554
 555   BOOST_CHECK_CLOSE(
 556   cdf(mypoisson, 2.),  // k events = 2.
 557     0.238103305553545, // probability.
 558       5e-13);
 559
 560
 561 #if 0
 562   // Compare cdf from finite sum of pdf and gamma_q.
 563   using boost::math::cdf;
 564   using boost::math::pdf;
 565
 566   double mean = 4.;
 567   cout.precision(17); // double max_digits10
 568   cout.setf(ios::showpoint);
 569   cout << showpoint << endl;  // Ensure trailing zeros are shown.
 570   // This also helps show the expected precision max_digits10
 571   //cout.unsetf(ios::showpoint); // No trailing zeros are shown.
 572
 573   cout << "k          pdf                     sum                  cdf                   diff" << endl;
 574   double sum = 0.;
 575   for (int i = 0; i <= 50; i++)
 576   {
 577    cout << i << ' ' ;
 578    double p =  pdf(poisson_distribution<double>(mean), static_cast<double>(i));
 579    sum += p;
 580
 581    cout << p << ' ' << sum << ' '
 582    << cdf(poisson_distribution<double>(mean), static_cast<double>(i)) << ' ';
 583      {
 584        cout << boost::math::gamma_q<double>(i+1, mean); // cdf
 585        double diff = boost::math::gamma_q<double>(i+1, mean) - sum; // cdf -sum
 586        cout << setprecision (2) << ' ' << diff; // 0 0 to 4, 1 eps 5 to 9, 10 to 20 2 eps, 21 upwards 3 eps
 587
 588      }
 589     BOOST_CHECK_CLOSE(
 590     cdf(mypoisson, static_cast<double>(i)),
 591       sum, // of pdfs.
 592       4e-14); // Fails at 2e-14
 593    // This call puts the precision etc back to default 6 !!!
 594    cout << setprecision(17) << showpoint;
 595
 596
 597      cout << endl;
 598   }
 599
 600    cout << cdf(poisson_distribution<double>(5), static_cast<double>(0)) << ' ' << endl; // 0.006737946999085467
 601    cout << cdf(poisson_distribution<double>(5), static_cast<double>(1)) << ' ' << endl; // 0.040427681994512805
 602    cout << cdf(poisson_distribution<double>(2), static_cast<double>(3)) << ' ' << endl; // 0.85712346049854715
 603
 604    { // Compare approximate formula in Wikipedia with quantile(half)
 605      for (int i = 1; i < 100; i++)
 606      {
 607        poisson_distribution<double> distn(static_cast<double>(i));
 608        cout << i << ' ' << median(distn) << ' ' << quantile(distn, 0.5) << ' '
 609          << median(distn) - quantile(distn, 0.5) << endl; // formula appears to be out-by-one??
 610      }  // so quantile(half) used via derived accressors.
 611    }
 612 #endif
 613
 614    // (Parameter value, arbitrarily zero, only communicates the floating-point type).
 615 #ifdef TEST_POISSON
 616   test_spots(0.0F); // Test float.
 617 #endif
 618 #ifdef TEST_DOUBLE
 619   test_spots(0.0); // Test double.
 620 #endif
 621 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 622   if (numeric_limits<long double>::digits10 > numeric_limits<double>::digits10)
 623   { // long double is better than double (so not MSVC where they are same).
 624 #ifdef TEST_LDOUBLE
 625      test_spots(0.0L); // Test long double.
 626 #endif
 627   }
 628
 629 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
 630 #ifdef TEST_REAL_CONCEPT
 631   test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
 632 #endif
 633 #endif
 634 #endif
 635
 636 } // BOOST_AUTO_TEST_CASE( test_main )
 637
 638 /*
 639
 640 Output:
 641
 642 Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_poisson.exe"
 643 Running 1 test case...
 644 mean(mypoisson, 4.) == 4.0000000000000000
 645 mean(mypoisson, 0.) == 4.0000000000000000
 646 cdf(mypoisson, 2.) == 0.23810330555354431
 647 pdf(mypoisson, 2.) == 0.14652511110987343
 648 *** No errors detected
 649
 650 */