test_poisson.cpp

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

     cdf(complement(poisson_distribution<RealType>(static_cast<RealType>(4.)), // mean
      static_cast<RealType>(0))),  // Zero k events (uses special case formula, not gamma).
      static_cast<RealType>(0.98168436111126578), // probability.
         tolerance);
  BOOST_CHECK_CLOSE( // Complement CDF
     cdf(complement(poisson_distribution<RealType>(static_cast<RealType>(1.)), // mean
      static_cast<RealType>(0))),  // Zero k events (uses special case formula, not gamma).
      static_cast<RealType>(0.63212055882855767), // probability.
         tolerance);

  // Example where k is bigger than max_factorial (>34 for float)
  // (therefore using log gamma so perhaps less accurate).
  BOOST_CHECK_CLOSE(
     cdf(poisson_distribution<RealType>(static_cast<RealType>(40.)), // mean
      static_cast<RealType>(40)),  // k events. 
      static_cast<RealType>(0.5419181783625430), // probability.
         tolerance);

   // Quantile & complement.
  BOOST_CHECK_CLOSE(
    boost::math::quantile(
         poisson_distribution<RealType>(5),  // mean.
         static_cast<RealType>(0.615960654833065)),  //  probability.
         static_cast<RealType>(5.), // Expect k = 5
         tolerance/5); // 

  // EQUAL is too optimistic - fails [5.0000000000000124 != 5]
  // BOOST_CHECK_EQUAL(boost::math::quantile( // 
  //       poisson_distribution<RealType>(5.),  // mean.
  //       static_cast<RealType>(0.615960654833065)),  //  probability.
  //       static_cast<RealType>(5.)); // Expect k = 5 events.
 
  BOOST_CHECK_CLOSE(boost::math::quantile(
         poisson_distribution<RealType>(4),  // mean.
         static_cast<RealType>(0.785130387030406)),  //  probability.
         static_cast<RealType>(5.), // Expect k = 5 events.
         tolerance/5); 

  // Check on quantile of other examples of inverse of cdf.
  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(boost::math::quantile( // inverse of cdf above.
         poisson_distribution<RealType>(10.),  // mean.
         static_cast<RealType>(0.5830397501929856)),  //  probability.
         static_cast<RealType>(10.), // Expect k = 10 events.
         tolerance/5); 


  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);

  BOOST_CHECK_CLOSE(boost::math::quantile( // inverse of cdf above.
         poisson_distribution<RealType>(4.),  // mean.
         static_cast<RealType>(0.785130387030406)),  //  probability.
         static_cast<RealType>(5.), // Expect k = 10 events.
         tolerance/5); 



  //BOOST_CHECK_CLOSE(boost::math::quantile(
  //       poisson_distribution<RealType>(5),  // mean.
  //       static_cast<RealType>(0.785130387030406)),  //  probability.
  //        // 6.1882832344329559 result but MathCAD givest smallest integer ppois(k, mean) >= prob
  //       static_cast<RealType>(6.), // Expect k = 6 events. 
  //       tolerance/5); 

  //BOOST_CHECK_CLOSE(boost::math::quantile(
  //       poisson_distribution<RealType>(5),  // mean.
  //       static_cast<RealType>(0.77)),  //  probability.
  //        // 6.1882832344329559 result but MathCAD givest smallest integer ppois(k, mean) >= prob
  //       static_cast<RealType>(7.), // Expect k = 6 events. 
  //       tolerance/5); 

  //BOOST_CHECK_CLOSE(boost::math::quantile(
  //       poisson_distribution<RealType>(5),  // mean.
  //       static_cast<RealType>(0.75)),  //  probability.
  //        // 6.1882832344329559 result but MathCAD givest smallest integer ppois(k, mean) >= prob
  //       static_cast<RealType>(6.), // Expect k = 6 events. 
  //       tolerance/5); 

  BOOST_CHECK_CLOSE(
    boost::math::quantile(
         complement(
           poisson_distribution<RealType>(4),
           static_cast<RealType>(1 - 0.785130387030406))),  // complement.
           static_cast<RealType>(5), // Expect k = 5 events.
         tolerance/5);

  BOOST_CHECK_EQUAL(boost::math::quantile( // Check case when probability < cdf(0) (== pdf(0))
         poisson_distribution<RealType>(1),  // mean is small, so cdf and pdf(0) are about 0.35.
         static_cast<RealType>(0.0001)),  //  probability < cdf(0).
         static_cast<RealType>(0)); // Expect k = 0 events exactly.
          
  BOOST_CHECK_EQUAL(
    boost::math::quantile(
         complement(
           poisson_distribution<RealType>(1),
           static_cast<RealType>(0.9999))),  // complement, so 1-probability < cdf(0)
           static_cast<RealType>(0)); // Expect k = 0 events exactly.

  //
  // Test quantile policies against test data:
  //
#define T RealType
#include "poisson_quantile.ipp"

  for(unsigned i = 0; i < poisson_quantile_data.size(); ++i)
  {
     using namespace boost::math::policies;
     typedef policy<discrete_quantile<real> > P1;
     typedef policy<discrete_quantile<integer_round_down> > P2;
     typedef policy<discrete_quantile<integer_round_up> > P3;
     typedef policy<discrete_quantile<integer_round_outwards> > P4;
     typedef policy<discrete_quantile<integer_round_inwards> > P5;
     typedef policy<discrete_quantile<integer_round_nearest> > P6;
     RealType tol = boost::math::tools::epsilon<RealType>() * 20;
     if(!boost::is_floating_point<RealType>::value)
        tol *= 7;
     //
     // Check full real value first:
     //
     poisson_distribution<RealType, P1> p1(poisson_quantile_data[i][0]);
     RealType x = quantile(p1, poisson_quantile_data[i][1]);
     BOOST_CHECK_CLOSE_FRACTION(x, poisson_quantile_data[i][2], tol);
     x = quantile(complement(p1, poisson_quantile_data[i][1]));
     BOOST_CHECK_CLOSE_FRACTION(x, poisson_quantile_data[i][3], tol);
     //
     // Now with round down to integer:
     //
     poisson_distribution<RealType, P2> p2(poisson_quantile_data[i][0]);
     x = quantile(p2, poisson_quantile_data[i][1]);
     BOOST_CHECK_EQUAL(x, floor(poisson_quantile_data[i][2]));
     x = quantile(complement(p2, poisson_quantile_data[i][1]));
     BOOST_CHECK_EQUAL(x, floor(poisson_quantile_data[i][3]));
     //
     // Now with round up to integer:
     //
     poisson_distribution<RealType, P3> p3(poisson_quantile_data[i][0]);
     x = quantile(p3, poisson_quantile_data[i][1]);
     BOOST_CHECK_EQUAL(x, ceil(poisson_quantile_data[i][2]));
     x = quantile(complement(p3, poisson_quantile_data[i][1]));
     BOOST_CHECK_EQUAL(x, ceil(poisson_quantile_data[i][3]));
     //
     // Now with round to integer "outside":
     //
     poisson_distribution<RealType, P4> p4(poisson_quantile_data[i][0]);
     x = quantile(p4, poisson_quantile_data[i][1]);
     BOOST_CHECK_EQUAL(x, poisson_quantile_data[i][1] < 0.5f ? floor(poisson_quantile_data[i][2]) : ceil(poisson_quantile_data[i][2]));
     x = quantile(complement(p4, poisson_quantile_data[i][1]));
     BOOST_CHECK_EQUAL(x, poisson_quantile_data[i][1] < 0.5f ? ceil(poisson_quantile_data[i][3]) : floor(poisson_quantile_data[i][3]));
     //
     // Now with round to integer "inside":
     //
     poisson_distribution<RealType, P5> p5(poisson_quantile_data[i][0]);
     x = quantile(p5, poisson_quantile_data[i][1]);
     BOOST_CHECK_EQUAL(x, poisson_quantile_data[i][1] < 0.5f ? ceil(poisson_quantile_data[i][2]) : floor(poisson_quantile_data[i][2]));
     x = quantile(complement(p5, poisson_quantile_data[i][1]));
     BOOST_CHECK_EQUAL(x, poisson_quantile_data[i][1] < 0.5f ? floor(poisson_quantile_data[i][3]) : ceil(poisson_quantile_data[i][3]));
     //
     // Now with round to nearest integer:
     //
     poisson_distribution<RealType, P6> p6(poisson_quantile_data[i][0]);
     x = quantile(p6, poisson_quantile_data[i][1]);
     BOOST_CHECK_EQUAL(x, floor(poisson_quantile_data[i][2] + 0.5f));
     x = quantile(complement(p6, poisson_quantile_data[i][1]));
     BOOST_CHECK_EQUAL(x, floor(poisson_quantile_data[i][3] + 0.5f));
  }

} // template <class RealType>void test_spots(RealType)

//

int test_main(int, char* [])
{
  // Check that can construct normal distribution using the two convenience methods:
  using namespace boost::math;
  poisson myp1(2); // Using typedef
   poisson_distribution<> myp2(2); // Using default RealType double.

   // Basic sanity-check spot values.

  // Some plain double examples & tests:
  cout.precision(17); // double max_digits10
  cout.setf(ios::showpoint);
  
  poisson mypoisson(4.); // // mean = 4, default FP type is double.
  cout << "mean(mypoisson, 4.) == " << mean(mypoisson) << endl;
  cout << "mean(mypoisson, 0.) == " << mean(mypoisson) << endl;
  cout << "cdf(mypoisson, 2.) == " << cdf(mypoisson, 2.) << endl;
  cout << "pdf(mypoisson, 2.) == " << pdf(mypoisson, 2.) << endl;
  
  // poisson mydudpoisson(0.);
  // throws (if BOOST_MATH_DOMAIN_ERROR_POLICY == throw_on_error).

 
  BOOST_CHECK_THROW(poisson mydudpoisson(-1), std::domain_error);// Mean must be > 0.
  BOOST_CHECK_THROW(poisson mydudpoisson(-1), std::logic_error);// Mean must be > 0.
  // Passes the check because logic_error is a parent????
  // BOOST_CHECK_THROW(poisson mydudpoisson(-1), std::overflow_error); // fails the check
  // because overflow_error is unrelated - except from std::exception
  BOOST_CHECK_THROW(cdf(mypoisson, -1), std::domain_error); // k must be >= 0

  BOOST_CHECK_EQUAL(mean(mypoisson), 4.);
  BOOST_CHECK_CLOSE(
  pdf(mypoisson, 2.),  // k events = 2. 
    1.465251111098740E-001, // probability.
      5e-13);

  BOOST_CHECK_CLOSE(
  cdf(mypoisson, 2.),  // k events = 2. 
    0.238103305553545, // probability.
      5e-13);


#if 0
  // Compare cdf from finite sum of pdf and gamma_q.
  using boost::math::cdf;
  using boost::math::pdf;

  double mean = 4.;
  cout.precision(17); // double max_digits10
  cout.setf(ios::showpoint);
  cout << showpoint << endl;  // Ensure trailing zeros are shown.
  // This also helps show the expected precision max_digits10
  //cout.unsetf(ios::showpoint); // No trailing zeros are shown.

  cout << "k          pdf                     sum                  cdf                   diff" << endl;
  double sum = 0.;
  for (int i = 0; i <= 50; i++)
  {
   cout << i << ' ' ;
   double p =  pdf(poisson_distribution<double>(mean), static_cast<double>(i));
   sum += p;

   cout << p << ' ' << sum << ' ' 
   << cdf(poisson_distribution<double>(mean), static_cast<double>(i)) << ' ';
     {
       cout << boost::math::gamma_q<double>(i+1, mean); // cdf
       double diff = boost::math::gamma_q<double>(i+1, mean) - sum; // cdf -sum
       cout << setprecision (2) << ' ' << diff; // 0 0 to 4, 1 eps 5 to 9, 10 to 20 2 eps, 21 upwards 3 eps
      
     }
    BOOST_CHECK_CLOSE(
    cdf(mypoisson, static_cast<double>(i)),
      sum, // of pdfs.
      4e-14); // Fails at 2e-14
   // This call puts the precision etc back to default 6 !!!
   cout << setprecision(17) << showpoint;


     cout << endl;
  }

   cout << cdf(poisson_distribution<double>(5), static_cast<double>(0)) << ' ' << endl; // 0.006737946999085467
   cout << cdf(poisson_distribution<double>(5), static_cast<double>(1)) << ' ' << endl; // 0.040427681994512805
   cout << cdf(poisson_distribution<double>(2), static_cast<double>(3)) << ' ' << endl; // 0.85712346049854715 

   { // Compare approximate formula in Wikipedia with quantile(half)
     for (int i = 1; i < 100; i++)
     {
       poisson_distribution<double> distn(static_cast<double>(i));
       cout << i << ' ' << median(distn) << ' ' << quantile(distn, 0.5) << ' ' 
         << median(distn) - quantile(distn, 0.5) << endl; // formula appears to be out-by-one??
     }  // so quantile(half) used via derived accressors.
   }
#endif

   // (Parameter value, arbitrarily zero, only communicates the floating-point type).
#ifdef TEST_POISSON
  test_spots(0.0F); // Test float.
#endif
#ifdef TEST_DOUBLE
  test_spots(0.0); // Test double.
#endif
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  if (numeric_limits<long double>::digits10 > numeric_limits<double>::digits10)
  { // long double is better than double (so not MSVC where they are same).
#ifdef TEST_LDOUBLE
     test_spots(0.0L); // Test long double.
#endif
  }

  #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
#ifdef TEST_REAL_CONCEPT
  test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
#endif
  #endif
#endif
   return 0;
} // int test_main(int, char* [])

/*

Output:

Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_poisson.exe"
Running 1 test case...
mean(mypoisson, 4.) == 4.0000000000000000
mean(mypoisson, 0.) == 4.0000000000000000
cdf(mypoisson, 2.) == 0.23810330555354431
pdf(mypoisson, 2.) == 0.14652511110987343
*** No errors detected

*/