test_negative_binomial.cpp

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

// test_negative_binomial.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)

// Tests for Negative Binomial Distribution.

// Note that these defines must be placed BEFORE #includes.
#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
// because several tests overflow & underflow by design.
#define BOOST_MATH_DISCRETE_QUANTILE_POLICY real

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

#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

#include <boost/math/concepts/real_concept.hpp> // for real_concept
using ::boost::math::concepts::real_concept;

#include <boost/math/distributions/negative_binomial.hpp> // for negative_binomial_distribution
using boost::math::negative_binomial_distribution;

#include <boost/math/special_functions/gamma.hpp>
  using boost::math::lgamma;  // log gamma

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

template <class RealType>
void test_spot( // Test a single spot value against 'known good' values.
               RealType N,    // Number of successes.
               RealType k,    // Number of failures.
               RealType p,    // Probability of success_fraction.
               RealType P,    // CDF probability.
               RealType Q,    // Complement of CDF.
               RealType tol)  // Test tolerance.
{
   boost::math::negative_binomial_distribution<RealType> bn(N, p);
   BOOST_CHECK_EQUAL(N, bn.successes());
   BOOST_CHECK_EQUAL(p, bn.success_fraction());
   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 guarantee that 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);
      }
    }
    // estimate success ratio:
    BOOST_CHECK_CLOSE(
      negative_binomial_distribution<RealType>::find_lower_bound_on_p(
      N+k, N, P),
      p, tol);
    // Note we bump up the sample size here, purely for the sake of the test,
    // internally the function has to adjust the sample size so that we get
    // the right upper bound, our test undoes this, so we can verify the result.
    BOOST_CHECK_CLOSE(
      negative_binomial_distribution<RealType>::find_upper_bound_on_p(
      N+k+1, N, Q),
      p, tol);

    if(Q < P)
    {
       //
       // We check two things here, that the upper and lower bounds
       // are the right way around, and that they do actually bracket
       // the naive estimate of p = successes / (sample size)
       //
      BOOST_CHECK(
        negative_binomial_distribution<RealType>::find_lower_bound_on_p(
        N+k, N, Q)
        <=
        negative_binomial_distribution<RealType>::find_upper_bound_on_p(
        N+k, N, Q)
        );
      BOOST_CHECK(
        negative_binomial_distribution<RealType>::find_lower_bound_on_p(
        N+k, N, Q)
        <=
        N / (N+k)
        );
      BOOST_CHECK(
        N / (N+k)
        <=
        negative_binomial_distribution<RealType>::find_upper_bound_on_p(
        N+k, N, Q)
        );
    }
    else
    {
       // As above but when P is small.
      BOOST_CHECK(
        negative_binomial_distribution<RealType>::find_lower_bound_on_p(
        N+k, N, P)
        <=
        negative_binomial_distribution<RealType>::find_upper_bound_on_p(
        N+k, N, P)
        );
      BOOST_CHECK(
        negative_binomial_distribution<RealType>::find_lower_bound_on_p(
        N+k, N, P)
        <=
        N / (N+k)
        );
      BOOST_CHECK(
        N / (N+k)
        <=
        negative_binomial_distribution<RealType>::find_upper_bound_on_p(
        N+k, N, P)
        );
    }

    // Estimate sample size:
    BOOST_CHECK_CLOSE(
      negative_binomial_distribution<RealType>::find_minimum_number_of_trials(
      k, p, P),
      N+k, tol);
    BOOST_CHECK_CLOSE(
      negative_binomial_distribution<RealType>::find_maximum_number_of_trials(
         k, p, Q),
      N+k, 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);

    // Complement is not possible since sum is to infinity.
  } //
} // test_spot

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 1000 eps expressed as a percent, or
  // 1000 eps of type double expressed as a percent, whichever
  // is the larger.

  RealType tolerance = (std::max)
    (boost::math::tools::epsilon<RealType>(),
    static_cast<RealType>(std::numeric_limits<double>::epsilon()));
  tolerance *= 100 * 1000;

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

  RealType tol1eps = boost::math::tools::epsilon<RealType>() * 2; // Very tight, suit exact values.
  //RealType tol2eps = boost::math::tools::epsilon<RealType>() * 2; // Tight, suit exact values.
  RealType tol5eps = boost::math::tools::epsilon<RealType>() * 5; // Wider 5 epsilon.
  cout << "Tolerance 5 eps = " << tol5eps << "%." << endl;

  // Sources of spot test values:

  // MathCAD defines pbinom(k, r, p) (at about 64-bit double precision, about 16 decimal digits)
  // returns pr(X , k) when random variable X has the binomial distribution with parameters r and p.
  // 0 <= k
  // r > 0
  // 0 <= p <= 1
  // P = pbinom(30, 500, 0.05) = 0.869147702104609

  // And functions.wolfram.com

  using boost::math::negative_binomial_distribution;
  using  ::boost::math::negative_binomial;
  using  ::boost::math::cdf;
  using  ::boost::math::pdf;

  // Test negative binomial using cdf spot values from MathCAD cdf = pnbinom(k, r, p).
  // These test quantiles and complements as well.

  test_spot(  // pnbinom(1,2,0.5) = 0.5
  static_cast<RealType>(2),   // successes r
  static_cast<RealType>(1),   // Number of failures, k
  static_cast<RealType>(0.5), // Probability of success as fraction, p
  static_cast<RealType>(0.5), // Probability of result (CDF), P
  static_cast<RealType>(0.5),  // complement CCDF Q = 1 - P
  tolerance);

  test_spot( // pbinom(0, 2, 0.25)
  static_cast<RealType>(2),    // successes r
  static_cast<RealType>(0),    // Number of failures, k
  static_cast<RealType>(0.25),
  static_cast<RealType>(0.0625),                    // Probability of result (CDF), P
  static_cast<RealType>(0.9375),                    // Q = 1 - P
  tolerance);

  test_spot(  // pbinom(48,8,0.25)
  static_cast<RealType>(8),     // successes r
  static_cast<RealType>(48),    // Number of failures, k
  static_cast<RealType>(0.25),                    // Probability of success, p
  static_cast<RealType>(9.826582228110670E-1),     // Probability of result (CDF), P
  static_cast<RealType>(1 - 9.826582228110670E-1),   // Q = 1 - P
  tolerance);

  test_spot(  // pbinom(2,5,0.4)
  static_cast<RealType>(5),     // successes r
  static_cast<RealType>(2),     // Number of failures, k
  static_cast<RealType>(0.4),                    // Probability of success, p
  static_cast<RealType>(9.625600000000020E-2),     // Probability of result (CDF), P
  static_cast<RealType>(1 - 9.625600000000020E-2),   // Q = 1 - P
  tolerance);

  test_spot(  // pbinom(10,100,0.9)
  static_cast<RealType>(100),     // successes r
  static_cast<RealType>(10),     // Number of failures, k
  static_cast<RealType>(0.9),                    // Probability of success, p
  static_cast<RealType>(4.535522887695670E-1),     // Probability of result (CDF), P
  static_cast<RealType>(1 - 4.535522887695670E-1),   // Q = 1 - P
  tolerance);

  test_spot(  // pbinom(1,100,0.991)
  static_cast<RealType>(100),     // successes r
  static_cast<RealType>(1),     // Number of failures, k
  static_cast<RealType>(0.991),                    // Probability of success, p
  static_cast<RealType>(7.693413044217000E-1),     // Probability of result (CDF), P
  static_cast<RealType>(1 - 7.693413044217000E-1),   // Q = 1 - P
  tolerance);

  test_spot(  // pbinom(10,100,0.991)
  static_cast<RealType>(100),     // successes r
  static_cast<RealType>(10),     // Number of failures, k
  static_cast<RealType>(0.991),                    // Probability of success, p
  static_cast<RealType>(9.999999940939000E-1),     // Probability of result (CDF), P
  static_cast<RealType>(1 - 9.999999940939000E-1),   // Q = 1 - P
  tolerance);

if(std::numeric_limits<RealType>::is_specialized)