test_poisson.cpp

Boost provides free peer-reviewed portable C++ source libraries. We emphasize libraries that work

// test_poisson.cpp

// Copyright Paul A. Bristow 2007.
// Copyright John Maddock 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)

// Basic sanity test for Poisson Cumulative Distribution Function.

#define BOOST_MATH_DISCRETE_QUANTILE_POLICY real

#if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
#  define TEST_FLOAT
#  define TEST_DOUBLE
#  define TEST_LDOUBLE
#  define TEST_REAL_CONCEPT
#endif

#ifdef _MSC_VER
#  pragma warning(disable: 4127) // conditional expression is constant.
#endif

#include <boost/test/included/test_exec_monitor.hpp> // Boost.Test
#include <boost/test/floating_point_comparison.hpp>

#include <boost/math/concepts/real_concept.hpp> // for real_concept
#include <boost/math/distributions/poisson.hpp>
    using boost::math::poisson_distribution;
#include <boost/math/tools/test.hpp> // for real_concept

#include <boost/math/special_functions/gamma.hpp> // for (incomplete) gamma.
//   using boost::math::qamma_Q;

#include <iostream>
   using std::cout;
   using std::endl;
   using std::setprecision;
   using std::showpoint;
   using std::ios;
#include <limits>
  using std::numeric_limits;

template <class RealType> // Any floating-point type RealType.
void test_spots(RealType)
{
  // Basic sanity checks, tolerance is about numeric_limits<RealType>::digits10 decimal places,
   // guaranteed for type RealType, eg 6 for float, 15 for double,
   // expressed as a percentage (so -2) for BOOST_CHECK_CLOSE,

   int decdigits = numeric_limits<RealType>::digits10;
  // May eb >15 for 80 and 128-bit FP typtes.
  if (decdigits <= 0)
  { // decdigits is not defined, for example real concept,
    // so assume precision of most test data is double (for example, MathCAD).
     decdigits = numeric_limits<double>::digits10; // == 15 for 64-bit
  }
  if (decdigits > 15 ) // numeric_limits<double>::digits10)
  { // 15 is the accuracy of the MathCAD test data.
    decdigits = 15; // numeric_limits<double>::digits10;
  }

   decdigits -= 1; // Perhaps allow some decimal digit(s) margin of numerical error.
   RealType tolerance = static_cast<RealType>(std::pow(10., static_cast<double>(2-decdigits))); // 1e-6 (-2 so as %)
   tolerance *= 2; // Allow some bit(s) small margin (2 means + or - 1 bit) of numerical error.
   // Typically 2e-13% = 2e-15 as fraction for double.

   // Sources of spot test values:

  // Many be some combinations for which the result is 'exact',
  // or at least is good to 40 decimal digits.
   // 40 decimal digits includes 128-bit significand User Defined Floating-Point types,
   
   // Best source of accurate values is:
   // Mathworld online calculator (40 decimal digits precision, suitable for up to 128-bit significands)
   // http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=GammaRegularized
   // GammaRegularized is same as gamma incomplete, gamma or gamma_q(a, x) or Q(a, z).

  // http://documents.wolfram.com/calculationcenter/v2/Functions/ListsMatrices/Statistics/PoissonDistribution.html

  // MathCAD defines ppois(k, lambda== mean) as k integer, k >=0.
  // ppois(0, 5) =  6.73794699908547e-3
  // ppois(1, 5) = 0.040427681994513;
  // ppois(10, 10) = 5.830397501929850E-001
  // ppois(10, 1) = 9.999999899522340E-001
  // ppois(5,5) = 0.615960654833065

  // qpois returns inverse Poission distribution, that is the smallest (floor) k so that ppois(k, lambda) >= p
  // p is real number, real mean lambda > 0
  // k is approximately the integer for which probability(X <= k) = p
  // when random variable X has the Poisson distribution with parameters lambda.
  // Uses discrete bisection.
  // qpois(6.73794699908547e-3, 5) = 1
  // qpois(0.040427681994513, 5) = 

  // Test Poisson with spot values from MathCAD 'known good'.

  using boost::math::poisson_distribution;
  using  ::boost::math::poisson;
  using  ::boost::math::cdf;
  using  ::boost::math::pdf;

   // Check that bad arguments throw.
   BOOST_CHECK_THROW(
   cdf(poisson_distribution<RealType>(static_cast<RealType>(0)), // mean zero is bad.
      static_cast<RealType>(0)),  // even for a good k.
      std::domain_error); // Expected error to be thrown.

    BOOST_CHECK_THROW(
   cdf(poisson_distribution<RealType>(static_cast<RealType>(-1)), // mean negative is bad.
      static_cast<RealType>(0)),
      std::domain_error);

   BOOST_CHECK_THROW(
   cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unit OK,
      static_cast<RealType>(-1)),  // but negative events is bad.
      std::domain_error);

  BOOST_CHECK_THROW(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(0)), // mean zero is bad.
      static_cast<RealType>(99999)),  // for any k events. 
      std::domain_error);
  
  BOOST_CHECK_THROW(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(0)), // mean zero is bad.
      static_cast<RealType>(99999)),  // for any k events. 
      std::domain_error);

  BOOST_CHECK_THROW(
     quantile(poisson_distribution<RealType>(static_cast<RealType>(0)), // mean zero.
      static_cast<RealType>(0.5)),  // probability OK. 
      std::domain_error);

  BOOST_CHECK_THROW(
     quantile(poisson_distribution<RealType>(static_cast<RealType>(-1)), 
      static_cast<RealType>(-1)),  // bad probability. 
      std::domain_error);

  BOOST_CHECK_THROW(
     quantile(poisson_distribution<RealType>(static_cast<RealType>(1)), 
      static_cast<RealType>(-1)),  // bad probability. 
      std::domain_error);

  // Check some test values.

  BOOST_CHECK_CLOSE( // mode
     mode(poisson_distribution<RealType>(static_cast<RealType>(4))), // mode = mean = 4.
      static_cast<RealType>(4), // mode.
         tolerance);

  //BOOST_CHECK_CLOSE( // mode
  //   median(poisson_distribution<RealType>(static_cast<RealType>(4))), // mode = mean = 4.
  //    static_cast<RealType>(4), // mode.
      //   tolerance);
  poisson_distribution<RealType> dist4(static_cast<RealType>(40));

  BOOST_CHECK_CLOSE( // median
     median(dist4), // mode = mean = 4. median = 40.328333333333333 
      quantile(dist4, static_cast<RealType>(0.5)), // 39.332839138842637
         tolerance);

  // PDF
  BOOST_CHECK_CLOSE(
     pdf(poisson_distribution<RealType>(static_cast<RealType>(4)), // mean 4.
      static_cast<RealType>(0)),   
      static_cast<RealType>(1.831563888873410E-002), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     pdf(poisson_distribution<RealType>(static_cast<RealType>(4)), // mean 4.
      static_cast<RealType>(2)),   
      static_cast<RealType>(1.465251111098740E-001), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     pdf(poisson_distribution<RealType>(static_cast<RealType>(20)), // mean big.
      static_cast<RealType>(1)),   //  k small
      static_cast<RealType>(4.122307244877130E-008), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     pdf(poisson_distribution<RealType>(static_cast<RealType>(4)), // mean 4.
      static_cast<RealType>(20)),   //  K>> mean 
      static_cast<RealType>(8.277463646553730E-009), // probability.
         tolerance);
  
  // CDF
  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
      static_cast<RealType>(0)),  // zero k events. 
      static_cast<RealType>(3.678794411714420E-1), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
      static_cast<RealType>(1)),  // one k event. 
      static_cast<RealType>(7.357588823428830E-1), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
      static_cast<RealType>(2)),  // two k events. 
      static_cast<RealType>(9.196986029286060E-1), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
      static_cast<RealType>(10)),  // two k events. 
      static_cast<RealType>(9.999999899522340E-1), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
      static_cast<RealType>(15)),  // two k events. 
      static_cast<RealType>(9.999999999999810E-1), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
      static_cast<RealType>(16)),  // two k events. 
      static_cast<RealType>(9.999999999999990E-1), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
      static_cast<RealType>(17)),  // two k events. 
      static_cast<RealType>(1.), // probability unity for double.
         tolerance);

  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
      static_cast<RealType>(33)),  // k events at limit for float unchecked_factorial table. 
      static_cast<RealType>(1.), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(100)), // mean 100.
      static_cast<RealType>(33)),  // k events at limit for float unchecked_factorial table. 
      static_cast<RealType>(6.328271240363390E-15), // probability is tiny.
         tolerance * static_cast<RealType>(2e11)); // 6.3495253382825722e-015 MathCAD
      // Note that there two tiny probability are much more different.

   BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(100)), // mean 100.
      static_cast<RealType>(34)),  // k events at limit for float unchecked_factorial table. 
      static_cast<RealType>(1.898481372109020E-14), // probability is tiny.
         tolerance*static_cast<RealType>(2e11)); //         1.8984813721090199e-014 MathCAD


 BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(33)), // mean = k
      static_cast<RealType>(33)),  // k events above limit for float unchecked_factorial table. 
      static_cast<RealType>(5.461191812386560E-1), // probability.
         tolerance);

 BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(33)), // mean = k-1
      static_cast<RealType>(34)),  // k events above limit for float unchecked_factorial table. 
      static_cast<RealType>(6.133535681502950E-1), // probability.
         tolerance);

 BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(1)), // mean unity.
      static_cast<RealType>(34)),  // k events above limit for float unchecked_factorial table. 
      static_cast<RealType>(1.), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(5.)), // mean
      static_cast<RealType>(5)),  // k events. 
      static_cast<RealType>(0.615960654833065), // probability.
         tolerance);
  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(5.)), // mean
      static_cast<RealType>(1)),  // k events. 
      static_cast<RealType>(0.040427681994512805), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(5.)), // mean
      static_cast<RealType>(0)),  // k events (uses special case formula, not gamma). 
      static_cast<RealType>(0.006737946999085467), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(1.)), // mean
      static_cast<RealType>(0)),  // k events (uses special case formula, not gamma). 
      static_cast<RealType>(0.36787944117144233), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(10.)), // mean
      static_cast<RealType>(10)),  // k events. 
      static_cast<RealType>(0.5830397501929856), // probability.
         tolerance);

  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(4.)), // mean
      static_cast<RealType>(5)),  // k events. 
      static_cast<RealType>(0.785130387030406), // probability.
         tolerance);

  // complement CDF
  BOOST_CHECK_CLOSE( // Complement CDF
     cdf(complement(poisson_distribution<RealType>(static_cast<RealType>(4.)), // mean
      static_cast<RealType>(5))),  // k events. 
      static_cast<RealType>(1 - 0.785130387030406), // probability.
         tolerance);

  BOOST_CHECK_CLOSE( // Complement CDF