test_binomial.cpp

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

// test_binomial.cpp

// Copyright John Maddock 2006.
// Copyright  Paul A. Bristow 2007.

// 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 Binomial 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/math/concepts/real_concept.hpp> // for real_concept
using ::boost::math::concepts::real_concept;

#include <boost/math/distributions/binomial.hpp> // for binomial_distribution
using boost::math::binomial_distribution;

#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;
#include <limits>
using std::numeric_limits;

template <class RealType>
void test_spot(
     RealType N,    // Number of trials
     RealType k,    // Number of successes
     RealType p,    // Probability of success
     RealType P,    // CDF
     RealType Q,    // Complement of CDF
     RealType tol)  // Test tolerance
{
   boost::math::binomial_distribution<RealType> bn(N, p);
   BOOST_CHECK_CLOSE(
      cdf(bn, k), P, tol);
   if((P < 0.99) && (Q < 0.99))
   {
      //
      // We can only check this if P is not too close to 1,
      // so that we can guarentee Q is free of error:
      //
      BOOST_CHECK_CLOSE(
         cdf(complement(bn, k)), Q, tol);
      if(k != 0)
      {
         BOOST_CHECK_CLOSE(
            quantile(bn, P), k, tol);
      }
      else
      {
         // Just check quantile is very small:
         if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent) && (boost::is_floating_point<RealType>::value))
         {
            // Limit where this is checked: if exponent range is very large we may
            // run out of iterations in our root finding algorithm.
            BOOST_CHECK(quantile(bn, P) < boost::math::tools::epsilon<RealType>() * 10);
         }
      }
      if(k != 0)
      {
         BOOST_CHECK_CLOSE(
            quantile(complement(bn, Q)), k, tol);
      }
      else
      {
         // Just check quantile is very small:
         if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent) && (boost::is_floating_point<RealType>::value))
         {
            // Limit where this is checked: if exponent range is very large we may
            // run out of iterations in our root finding algorithm.
            BOOST_CHECK(quantile(complement(bn, Q)) < boost::math::tools::epsilon<RealType>() * 10);
         }
      }
      if(k > 0)
      {
         // estimate success ratio:
         // Note lower bound uses a different formual internally
         // from upper bound, have to adjust things to prevent
         // fencepost errors:
         BOOST_CHECK_CLOSE(
            binomial_distribution<RealType>::find_lower_bound_on_p(
               N, k+1, Q),
            p, tol);
         BOOST_CHECK_CLOSE(
            binomial_distribution<RealType>::find_upper_bound_on_p(
               N, k, P),
            p, tol);

         if(Q < P)
         {
            // Default method (Clopper Pearson)
            BOOST_CHECK(
               binomial_distribution<RealType>::find_lower_bound_on_p(
                  N, k, Q)
                  <=
               binomial_distribution<RealType>::find_upper_bound_on_p(
                  N, k, Q)
                  );
            BOOST_CHECK((
               binomial_distribution<RealType>::find_lower_bound_on_p(
                  N, k, Q)
                  <= k/N) && (k/N <=
               binomial_distribution<RealType>::find_upper_bound_on_p(
                  N, k, Q))
                  );
            // Bayes Method (Jeffreys Prior)
            BOOST_CHECK(
               binomial_distribution<RealType>::find_lower_bound_on_p(
               N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval)
                  <=
               binomial_distribution<RealType>::find_upper_bound_on_p(
                  N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval)
                  );
            BOOST_CHECK((
               binomial_distribution<RealType>::find_lower_bound_on_p(
                  N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval)
                  <= k/N) && (k/N <=
               binomial_distribution<RealType>::find_upper_bound_on_p(
                  N, k, Q, binomial_distribution<RealType>::jeffreys_prior_interval))
                  );
         }
         else
         {
            // Default method (Clopper Pearson)
            BOOST_CHECK(
               binomial_distribution<RealType>::find_lower_bound_on_p(
                  N, k, P)
                  <=
               binomial_distribution<RealType>::find_upper_bound_on_p(
                  N, k, P)
                  );
            BOOST_CHECK(
               (binomial_distribution<RealType>::find_lower_bound_on_p(
                  N, k, P)
                  <= k / N) && (k/N <=
               binomial_distribution<RealType>::find_upper_bound_on_p(
                  N, k, P))
                  );
            // Bayes Method (Jeffreys Prior)
            BOOST_CHECK(
               binomial_distribution<RealType>::find_lower_bound_on_p(
                  N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval)
                  <=
               binomial_distribution<RealType>::find_upper_bound_on_p(
                  N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval)
                  );
            BOOST_CHECK(
               (binomial_distribution<RealType>::find_lower_bound_on_p(
                  N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval)
                  <= k / N) && (k/N <=
               binomial_distribution<RealType>::find_upper_bound_on_p(
                  N, k, P, binomial_distribution<RealType>::jeffreys_prior_interval))
                  );
         }
      }
      //
      // estimate sample size:
      //
      BOOST_CHECK_CLOSE(
         binomial_distribution<RealType>::find_minimum_number_of_trials(
            k, p, P),
         N, tol);
      BOOST_CHECK_CLOSE(
         binomial_distribution<RealType>::find_maximum_number_of_trials(
            k, p, Q),
         N, tol);
   }

   // Double check consistency of CDF and PDF by computing
   // the finite sum:
   RealType sum = 0;
   for(unsigned i = 0; i <= k; ++i)
      sum += pdf(bn, RealType(i));
   BOOST_CHECK_CLOSE(
      sum, P, tol);
   // And complement as well:
   sum = 0;
   for(RealType i = N; i > k; i -= 1)
      sum += pdf(bn, i);
   if(P < 0.99)
   {
      BOOST_CHECK_CLOSE(
         sum, Q, tol);
   }
   else
   {
      // Not enough information content in P for Q to be meaningful
      RealType tol = (std::max)(2 * Q, boost::math::tools::epsilon<RealType>());
      BOOST_CHECK(sum < tol);
   }
}

template <class RealType> // Any floating-point type RealType.
void test_spots(RealType)
{
  // Basic sanity checks, test data is to double precision only
  // so set tolerance to 100eps expressed as a persent, or
  // 100eps of type double expressed as a persent, whichever
  // is the larger.

  RealType tolerance = (std::max)
      (boost::math::tools::epsilon<RealType>(),
      static_cast<RealType>(std::numeric_limits<double>::epsilon()));
  tolerance *= 100 * 1000;
  RealType tol2 = boost::math::tools::epsilon<RealType>() * 5 * 100;  // 5 eps as a persent

  cout << "Tolerance = " << tolerance << "%." << endl;

  // Sources of spot test values:

  // MathCAD defines pbinom(k, n, p)
  // returns pr(X ,=k) when random variable X has the binomial distribution with parameters n and p.
  // 0 <= k ,= n
  // 0 <= p <= 1
  // P = pbinom(30, 500, 0.05) = 0.869147702104609

  using boost::math::binomial_distribution;
  using  ::boost::math::cdf;
  using  ::boost::math::pdf;

#if !defined(TEST_ROUNDING) || (TEST_ROUNDING == 0)
  // Test binomial using cdf spot values from MathCAD.
  // These test quantiles and complements as well.
  test_spot(
     static_cast<RealType>(500),                     // Sample size, N
     static_cast<RealType>(30),                      // Number of successes, k
     static_cast<RealType>(0.05),                    // Probability of success, p
     static_cast<RealType>(0.869147702104609),       // Probability of result (CDF), P
     static_cast<RealType>(1 - 0.869147702104609),   // Q = 1 - P
     tolerance);

  test_spot(
     static_cast<RealType>(500),                     // Sample size, N
     static_cast<RealType>(250),                     // Number of successes, k
     static_cast<RealType>(0.05),                    // Probability of success, p
     static_cast<RealType>(1),                       // Probability of result (CDF), P
     static_cast<RealType>(0),   // Q = 1 - P
     tolerance);