test_negative_binomial.cpp

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

{ // An extreme value test that takes 3 minutes using the real concept type
  // for which numeric_limits<RealType>::is_specialized == false, deliberately
  // and for which there is no Lanczos approximation defined (also deliberately)
  // giving a very slow computation, but with acceptable accuracy.
  // A possible enhancement might be to use a normal approximation for
  // extreme values, but this is not implemented.
  test_spot(  // pbinom(100000,100,0.001)
  static_cast<RealType>(100),     // successes r
  static_cast<RealType>(100000),     // Number of failures, k
  static_cast<RealType>(0.001),                    // Probability of success, p
  static_cast<RealType>(5.173047534260320E-1),     // Probability of result (CDF), P
  static_cast<RealType>(1 - 5.173047534260320E-1),   // Q = 1 - P
  tolerance*1000); // *1000 is OK 0.51730475350664229  versus

  // functions.wolfram.com
  //   for I[0.001](100, 100000+1) gives:
  // Wolfram       0.517304753506834882009032744488738352004003696396461766326713
  // JM nonLanczos 0.51730475350664229 differs at the 13th decimal digit.
  // MathCAD       0.51730475342603199 differs at 10th decimal digit.
}
 // End of single spot tests using RealType


  // Tests on PDF:
  BOOST_CHECK_CLOSE(
  pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(0.5)),
  static_cast<RealType>(0) ),  // k = 0.
  static_cast<RealType>(0.25), // 0
  tolerance);

  BOOST_CHECK_CLOSE(
  pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(4), static_cast<RealType>(0.5)),
  static_cast<RealType>(0)),  // k = 0.
  static_cast<RealType>(0.0625), // exact 1/16
  tolerance);

  BOOST_CHECK_CLOSE(
  pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
  static_cast<RealType>(0)),  // k = 0
  static_cast<RealType>(9.094947017729270E-13), // pbinom(0,20,0.25) = 9.094947017729270E-13
  tolerance);

  BOOST_CHECK_CLOSE(
  pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.2)),
  static_cast<RealType>(0)),  // k = 0
  static_cast<RealType>(1.0485760000000003e-014), // MathCAD 1.048576000000000E-14
  tolerance);

  BOOST_CHECK_CLOSE(
  pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(10), static_cast<RealType>(0.1)),
  static_cast<RealType>(0)),  // k = 0.
  static_cast<RealType>(1e-10), // MathCAD says zero, but suffers cancellation error?
  tolerance);

  BOOST_CHECK_CLOSE(
  pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.1)),
  static_cast<RealType>(0)),  // k = 0.
  static_cast<RealType>(1e-20), // MathCAD says zero, but suffers cancellation error?
  tolerance);


  BOOST_CHECK_CLOSE( // .
  pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.9)),
  static_cast<RealType>(0)),  // k.
  static_cast<RealType>(1.215766545905690E-1), // k=20  p = 0.9
  tolerance);

  // Tests on cdf:
  // MathCAD pbinom k, r, p) == failures, successes, probability.

  BOOST_CHECK_CLOSE(cdf(
    negative_binomial_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(0.5)), // successes = 2,prob 0.25
    static_cast<RealType>(0) ), // k = 0
    static_cast<RealType>(0.25), // probability 1/4
    tolerance);

  BOOST_CHECK_CLOSE(cdf(complement(
    negative_binomial_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(0.5)), // successes = 2,prob 0.25
    static_cast<RealType>(0) )), // k = 0
    static_cast<RealType>(0.75), // probability 3/4
    tolerance);
  BOOST_CHECK_CLOSE( // k = 1.
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
  static_cast<RealType>(1)),  // k =1.
  static_cast<RealType>(1.455191522836700E-11),
  tolerance);

  BOOST_CHECK_SMALL( // Check within an epsilon with CHECK_SMALL
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
  static_cast<RealType>(1)) -
  static_cast<RealType>(1.455191522836700E-11),
  tolerance );

  // Some exact (probably - judging by trailing zeros) values.
  BOOST_CHECK_CLOSE(
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
  static_cast<RealType>(0)),  // k.
  static_cast<RealType>(1.525878906250000E-5),
  tolerance);

  BOOST_CHECK_CLOSE(
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
  static_cast<RealType>(0)),  // k.
  static_cast<RealType>(1.525878906250000E-5),
  tolerance);

  BOOST_CHECK_SMALL(
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
  static_cast<RealType>(0)) -
  static_cast<RealType>(1.525878906250000E-5),
  tolerance );

  BOOST_CHECK_CLOSE( // k = 1.
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
  static_cast<RealType>(1)),  // k.
  static_cast<RealType>(1.068115234375010E-4),
  tolerance);

  BOOST_CHECK_CLOSE( // k = 2.
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
  static_cast<RealType>(2)),  // k.
  static_cast<RealType>(4.158020019531300E-4),
  tolerance);

  BOOST_CHECK_CLOSE( // k = 3.
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
  static_cast<RealType>(3)),  // k.bristow
  static_cast<RealType>(1.188278198242200E-3),
  tolerance);

  BOOST_CHECK_CLOSE( // k = 4.
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
  static_cast<RealType>(4)),  // k.
  static_cast<RealType>(2.781510353088410E-3),
  tolerance);

  BOOST_CHECK_CLOSE( // k = 5.
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
  static_cast<RealType>(5)),  // k.
  static_cast<RealType>(5.649328231811500E-3),
  tolerance);

  BOOST_CHECK_CLOSE( // k = 6.
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
  static_cast<RealType>(6)),  // k.
  static_cast<RealType>(1.030953228473680E-2),
  tolerance);

  BOOST_CHECK_CLOSE( // k = 7.
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
  static_cast<RealType>(7)),  // k.
  static_cast<RealType>(1.729983836412430E-2),
  tolerance);

  BOOST_CHECK_CLOSE( // k = 8.
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
  static_cast<RealType>(8)),  // k = n.
  static_cast<RealType>(2.712995628826370E-2),
  tolerance);

  BOOST_CHECK_CLOSE( //
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
  static_cast<RealType>(48)),  // k
  static_cast<RealType>(9.826582228110670E-1),
  tolerance);

  BOOST_CHECK_CLOSE( //
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
  static_cast<RealType>(64)),  // k
  static_cast<RealType>(9.990295004935590E-1),
  tolerance);

  BOOST_CHECK_CLOSE( //
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(5), static_cast<RealType>(0.4)),
  static_cast<RealType>(26)),  // k
  static_cast<RealType>(9.989686246611190E-1),
  tolerance);

  BOOST_CHECK_CLOSE( //
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(5), static_cast<RealType>(0.4)),
  static_cast<RealType>(2)),  // k failures
  static_cast<RealType>(9.625600000000020E-2),
  tolerance);

  BOOST_CHECK_CLOSE( //
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(50), static_cast<RealType>(0.9)),
  static_cast<RealType>(20)),  // k
  static_cast<RealType>(9.999970854144170E-1),
  tolerance);

  BOOST_CHECK_CLOSE( //
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(500), static_cast<RealType>(0.7)),
  static_cast<RealType>(200)),  // k
  static_cast<RealType>(2.172846379930550E-1),
  tolerance* 2);

  BOOST_CHECK_CLOSE( //
  cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(50), static_cast<RealType>(0.7)),
  static_cast<RealType>(20)),  // k
  static_cast<RealType>(4.550203671301790E-1),
  tolerance);

  // Tests of other functions, mean and other moments ...

  negative_binomial_distribution<RealType> dist(static_cast<RealType>(8), static_cast<RealType>(0.25));
  using namespace std; // ADL of std names.
  // mean:
  BOOST_CHECK_CLOSE(
    mean(dist), static_cast<RealType>(8 * (1 - 0.25) /0.25), tol5eps);
  BOOST_CHECK_CLOSE(
    mode(dist), static_cast<RealType>(21), tol1eps);
  // variance:
  BOOST_CHECK_CLOSE(
    variance(dist), static_cast<RealType>(8 * (1 - 0.25) / (0.25 * 0.25)), tol5eps);
  // std deviation:
  BOOST_CHECK_CLOSE(
    standard_deviation(dist), // 9.79795897113271239270
    static_cast<RealType>(9.797958971132712392789136298823565567864L), // using functions.wolfram.com
    //                              9.79795897113271152534  == sqrt(8 * (1 - 0.25) / (0.25 * 0.25)))
    tol5eps * 100);
  BOOST_CHECK_CLOSE(
    skewness(dist), //
    static_cast<RealType>(0.71443450831176036),
    // using http://mathworld.wolfram.com/skewness.html
    tolerance);
  BOOST_CHECK_CLOSE(
    kurtosis_excess(dist), //
    static_cast<RealType>(0.7604166666666666666666666666666666666666L), // using Wikipedia Kurtosis(excess) formula
    tol5eps * 100);
  BOOST_CHECK_CLOSE(
    kurtosis(dist), // true 
    static_cast<RealType>(3.76041666666666666666666666666666666666666L), // 
    tol5eps * 100);
  // hazard:
  RealType x = static_cast<RealType>(0.125);
  BOOST_CHECK_CLOSE(
  hazard(dist, x)
  , pdf(dist, x) / cdf(complement(dist, x)), tol5eps);
  // cumulative hazard:
  BOOST_CHECK_CLOSE(
  chf(dist, x), -log(cdf(complement(dist, x))), tol5eps);
  // coefficient_of_variation:
  BOOST_CHECK_CLOSE(
  coefficient_of_variation(dist)
  , standard_deviation(dist) / mean(dist), tol5eps);

  // Special cases for PDF:
  BOOST_CHECK_EQUAL(
  pdf(
  negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)), //
  static_cast<RealType>(0)),
  static_cast<RealType>(0) );

  BOOST_CHECK_EQUAL(
  pdf(
  negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)),
  static_cast<RealType>(0.0001)),
  static_cast<RealType>(0) );

  BOOST_CHECK_EQUAL(
  pdf(
  negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1)),
  static_cast<RealType>(0.001)),
  static_cast<RealType>(0) );

  BOOST_CHECK_EQUAL(
  pdf(
  negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1)),
  static_cast<RealType>(8)),
  static_cast<RealType>(0) );

  BOOST_CHECK_SMALL(
  pdf(
   negative_binomial_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(0.25)),
  static_cast<RealType>(0))-
  static_cast<RealType>(0.0625),
  2 * boost::math::tools::epsilon<RealType>() ); // Expect exact, but not quite.
  // numeric_limits<RealType>::epsilon()); // Not suitable for real concept!

  // Quantile boundary cases checks:
  BOOST_CHECK_EQUAL(
  quantile(  // zero P < cdf(0) so should be exactly zero.
  negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
  static_cast<RealType>(0)),