mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
706 lines
25 KiB
C++
706 lines
25 KiB
C++
// test_negative_binomial.cpp
|
|
|
|
// Copyright Paul A. Bristow 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.
|
|
|
|
#define BOOST_MATH_THROW_ON_DOMAIN_ERROR
|
|
|
|
#ifdef _MSC_VER
|
|
# pragma warning(disable: 4127) // conditional expression is constant.
|
|
# pragma warning(disable: 4100) // unreferenced formal parameter.
|
|
# pragma warning(disable: 4512) // assignment operator could not be generated.
|
|
#endif
|
|
|
|
#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/math/concepts/real_concept.hpp> // for real_concept
|
|
using ::boost::math::concepts::real_concept;
|
|
|
|
#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;
|
|
|
|
// Need all functions working to use this test_spot!
|
|
|
|
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.success());
|
|
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);
|
|
}
|
|
}
|
|
#if 0
|
|
// estimate success ratio:
|
|
BOOST_CHECK_CLOSE(
|
|
negative_binomial_distribution<RealType>::estimate_lower_bound_on_p(
|
|
N, k, Q),
|
|
p, tol);
|
|
BOOST_CHECK_CLOSE(
|
|
negative_binomial_distribution<RealType>::estimate_upper_bound_on_p(
|
|
N, k, P),
|
|
p, tol);
|
|
|
|
if(Q < P)
|
|
{
|
|
BOOST_CHECK(
|
|
negative_binomial_distribution<RealType>::estimate_lower_bound_on_p(
|
|
N, k, Q)
|
|
<
|
|
negative_binomial_distribution<RealType>::estimate_upper_bound_on_p(
|
|
N, k, Q)
|
|
);
|
|
}
|
|
else
|
|
{
|
|
BOOST_CHECK(
|
|
negative_binomial_distribution<RealType>::estimate_lower_bound_on_p(
|
|
N, k, P)
|
|
<
|
|
negative_binomial_distribution<RealType>::estimate_upper_bound_on_p(
|
|
N, k, P)
|
|
);
|
|
}
|
|
#endif
|
|
//
|
|
// Estimate sample size:
|
|
//
|
|
//BOOST_CHECK_CLOSE(
|
|
// negative_binomial_distribution<RealType>::estimate_number_of_successes(
|
|
// k, p, P),
|
|
// N, tol);
|
|
//BOOST_CHECK_CLOSE(
|
|
// negative_binomial_distribution<RealType>::estimate_number_of_successes(
|
|
// boost::math::complement(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);
|
|
}
|
|
} //
|
|
} // 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 100eps expressed as a percent, or
|
|
// 100eps 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;
|
|
|
|
// Sources of spot test values:
|
|
|
|
// MathCAD defines pbinom(k, r, p)
|
|
// 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
|
|
|
|
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.
|
|
// These test quantiles and complements as well.
|
|
//test_spot(
|
|
//static_cast<RealType>(2), // Trial n
|
|
//static_cast<RealType>(1), // Number of failures, k
|
|
//static_cast<RealType>(0.5), // Probability of success, p
|
|
//static_cast<RealType>(0.25), // Probability of result (CDF), P
|
|
//static_cast<RealType>(0.75), // 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);
|
|
|
|
//test_spot(
|
|
//static_cast<RealType>(500), // Sample size, N
|
|
//static_cast<RealType>(470), // Number of successes, k
|
|
//static_cast<RealType>(0.95), // Probability of success, p
|
|
//static_cast<RealType>(0.176470742656766), // Probability of result (CDF), P
|
|
//static_cast<RealType>(1 - 0.176470742656766), // Q = 1 - P
|
|
//tolerance * 10); // Note higher tolerance on this test!
|
|
|
|
//test_spot(
|
|
//static_cast<RealType>(500), // Sample size, N
|
|
//static_cast<RealType>(400), // 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);
|
|
|
|
//test_spot(
|
|
//static_cast<RealType>(500), // Sample size, N
|
|
//static_cast<RealType>(400), // Number of successes, k
|
|
//static_cast<RealType>(0.9), // Probability of success, p
|
|
//static_cast<RealType>(1.80180425681923E-11), // Probability of result (CDF), P
|
|
//static_cast<RealType>(1 - 1.80180425681923E-11), // Q = 1 - P
|
|
//tolerance);
|
|
|
|
//test_spot(
|
|
//static_cast<RealType>(500), // Sample size, N
|
|
//static_cast<RealType>(5), // Number of successes, k
|
|
//static_cast<RealType>(0.05), // Probability of success, p
|
|
//static_cast<RealType>(9.181808267643E-7), // Probability of result (CDF), P
|
|
//static_cast<RealType>(1 - 9.181808267643E-7), // Q = 1 - P
|
|
//tolerance);
|
|
|
|
//test_spot(
|
|
//static_cast<RealType>(2), // Sample size, N
|
|
//static_cast<RealType>(1), // Number of successes, k
|
|
//static_cast<RealType>(0.5), // Probability of success, p
|
|
//static_cast<RealType>(0.75), // Probability of result (CDF), P
|
|
//static_cast<RealType>(0.25), // Q = 1 - P
|
|
//tolerance);
|
|
|
|
//test_spot(
|
|
//static_cast<RealType>(8), // Sample size, N
|
|
//static_cast<RealType>(3), // Number of successes, k
|
|
//static_cast<RealType>(0.25), // Probability of success, p
|
|
//static_cast<RealType>(0.8861846923828125), // Probability of result (CDF), P
|
|
//static_cast<RealType>(1 - 0.8861846923828125), // Q = 1 - P
|
|
//tolerance);
|
|
|
|
//test_spot(
|
|
//static_cast<RealType>(8), // Sample size, N
|
|
//static_cast<RealType>(0), // Number of successes, k
|
|
//static_cast<RealType>(0.25), // Probability of success, p
|
|
//static_cast<RealType>(0.1001129150390625), // Probability of result (CDF), P
|
|
//static_cast<RealType>(1 - 0.1001129150390625), // Q = 1 - P
|
|
//tolerance);
|
|
|
|
//test_spot(
|
|
//static_cast<RealType>(8), // Sample size, N
|
|
//static_cast<RealType>(1), // Number of successes, k
|
|
//static_cast<RealType>(0.25), // Probability of success, p
|
|
//static_cast<RealType>(0.36708068847656244), // Probability of result (CDF), P
|
|
//static_cast<RealType>(1 - 0.36708068847656244), // Q = 1 - P
|
|
//tolerance);
|
|
|
|
//test_spot(
|
|
//static_cast<RealType>(8), // Sample size, N
|
|
//static_cast<RealType>(4), // Number of successes, k
|
|
//static_cast<RealType>(0.25), // Probability of success, p
|
|
//static_cast<RealType>(0.9727020263671875), // Probability of result (CDF), P
|
|
//static_cast<RealType>(1 - 0.9727020263671875), // Q = 1 - P
|
|
//tolerance);
|
|
|
|
//test_spot(
|
|
//static_cast<RealType>(8), // Sample size, N
|
|
//static_cast<RealType>(7), // Number of successes, k
|
|
//static_cast<RealType>(0.25), // Probability of success, p
|
|
//static_cast<RealType>(0.9999847412109375), // Probability of result (CDF), P
|
|
//static_cast<RealType>(1 - 0.9999847412109375), // Q = 1 - P
|
|
//tolerance);
|
|
|
|
// End of single spot tests using RealType
|
|
|
|
// Tests on cdf:
|
|
// MathCAD pbinom k, r, p) == failures, successes,
|
|
|
|
BOOST_CHECK_CLOSE(cdf( // pnbinom(1,2,0.5) = 0.5
|
|
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);
|
|
|
|
// 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.
|
|
static_cast<RealType>(0.25), // 0
|
|
tolerance);
|
|
|
|
|
|
BOOST_CHECK_CLOSE(
|
|
pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(0.5)),
|
|
static_cast<RealType>(1)), // k.
|
|
static_cast<RealType>(0.25), //
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE(
|
|
cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(10)), // k = 10
|
|
static_cast<RealType>(1.820643361830440E-6), // pbinom(10,20,0.25) = 1.820643361830440E-6
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE( // k = 0 use different formula - only exp so more accurate.
|
|
cdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(0)), // k.
|
|
static_cast<RealType>(9.094947017729270E-13), // k=0 p = 0.25
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE( // k = 20 use different formula - only exp so more accurate.
|
|
pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(20)), // k == n.
|
|
static_cast<RealType>(0.00000000000090949470177292823791), // k=20 p = 0.25
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE( // k = 1.
|
|
pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(20), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(1)), // k.
|
|
static_cast<RealType>(0.02114141292622662149369716644287), // k=1 p = 0.25
|
|
tolerance);
|
|
|
|
// Some exact (probably - judging by trailing zeros) values.
|
|
BOOST_CHECK_CLOSE(
|
|
pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(0)), // k.
|
|
static_cast<RealType>(0.10011291503906250000000000000000), // k=0 p = 0.25
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE( // k = 1.
|
|
pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(1)), // k.
|
|
static_cast<RealType>(0.26696777343750000000000000000000), // k=1 p = 0.25
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE( // k = 2.
|
|
pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(2)), // k.
|
|
static_cast<RealType>(0.31146240234375000000000000000000), // k=2 p = 0.25
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE( // k = 3.
|
|
pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(3)), // k.
|
|
static_cast<RealType>(0.20764160156250000000000000000000), // k=3 p = 0.25
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE( // k = 7.
|
|
pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(7)), // k.
|
|
static_cast<RealType>(0.00036621093750000000000000000000), // k=7 p = 0.25
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE( // k = 8.
|
|
pdf(negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(8)), // k = n.
|
|
static_cast<RealType>(0.00001525878906250000000000000000), // k=8 p = 0.25
|
|
tolerance);
|
|
|
|
RealType tol2 = boost::math::tools::epsilon<RealType>() * 5;
|
|
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), tol2);
|
|
// variance:
|
|
BOOST_CHECK_CLOSE(
|
|
variance(dist)
|
|
, static_cast<RealType>(8 * (1-0.25) / (0.25 * 0.25)), tol2);
|
|
// std deviation:
|
|
BOOST_CHECK_CLOSE(
|
|
standard_deviation(dist)
|
|
, static_cast<RealType>(sqrt(8 * (1-0.25) / (0.25 * 0.25))), tol2);
|
|
// hazard:
|
|
RealType x = static_cast<RealType>(0.125);
|
|
BOOST_CHECK_CLOSE(
|
|
hazard(dist, x)
|
|
, pdf(dist, x) / cdf(complement(dist, x)), tol2);
|
|
// cumulative hazard:
|
|
BOOST_CHECK_CLOSE(
|
|
chf(dist, x)
|
|
, -log(cdf(complement(dist, x))), tol2);
|
|
// coefficient_of_variation:
|
|
BOOST_CHECK_CLOSE(
|
|
coefficient_of_variation(dist)
|
|
, standard_deviation(dist) / mean(dist), tol2);
|
|
|
|
// 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) // Expected pdf.
|
|
);
|
|
BOOST_CHECK_EQUAL(
|
|
pdf(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)),
|
|
static_cast<RealType>(0.0001)), // Expected pdf.
|
|
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)), // Expected pdf.
|
|
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) // Expected pdf.
|
|
);
|
|
BOOST_CHECK_EQUAL(
|
|
pdf(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(0)),
|
|
static_cast<RealType>(1) // Expected pdf.
|
|
);
|
|
|
|
// Check that duff arguments throw domain_error:
|
|
BOOST_CHECK_THROW(
|
|
pdf(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(-1), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(0)), std::domain_error
|
|
);
|
|
BOOST_CHECK_THROW(
|
|
pdf(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
|
|
static_cast<RealType>(0)), std::domain_error
|
|
);
|
|
BOOST_CHECK_THROW(
|
|
pdf(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
|
|
static_cast<RealType>(0)), std::domain_error
|
|
);
|
|
BOOST_CHECK_THROW(
|
|
pdf(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(-1)), std::domain_error
|
|
);
|
|
BOOST_CHECK_THROW(
|
|
pdf(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(9)), std::domain_error
|
|
);
|
|
BOOST_CHECK_THROW(
|
|
cdf(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(-1)), std::domain_error
|
|
);
|
|
BOOST_CHECK_THROW(
|
|
cdf(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(9)), std::domain_error
|
|
);
|
|
BOOST_CHECK_THROW(
|
|
cdf(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
|
|
static_cast<RealType>(0)), std::domain_error
|
|
);
|
|
BOOST_CHECK_THROW(
|
|
cdf(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
|
|
static_cast<RealType>(0)), std::domain_error
|
|
);
|
|
BOOST_CHECK_THROW(
|
|
quantile(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
|
|
static_cast<RealType>(0)), std::domain_error
|
|
);
|
|
BOOST_CHECK_THROW(
|
|
quantile(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
|
|
static_cast<RealType>(0)), std::domain_error
|
|
);
|
|
BOOST_CHECK_EQUAL(
|
|
cdf(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
|
|
static_cast<RealType>(8)), static_cast<RealType>(1)
|
|
);
|
|
BOOST_CHECK_EQUAL(
|
|
cdf(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0)),
|
|
static_cast<RealType>(7)), static_cast<RealType>(1)
|
|
);
|
|
BOOST_CHECK_EQUAL(
|
|
cdf(
|
|
negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1)),
|
|
static_cast<RealType>(7)), static_cast<RealType>(0)
|
|
);
|
|
// End of single checks on exact and check throwing 'duff' out-of-domain values.
|
|
} // template <class RealType> void test_spots(RealType) // Any floating-point type RealType.
|
|
|
|
int test_main(int, char* [])
|
|
{
|
|
// Check that can generate negative_binomial distribution using the two convenience methods:
|
|
using namespace boost::math;
|
|
negative_binomial mynb1(2., 0.5); // Using typedef - default type is double.
|
|
negative_binomial_distribution<> myf2(2., 0.5); // Using default RealType double.
|
|
|
|
// Basic sanity-check spot values.
|
|
#ifdef BOOST_MATH_THROW_ON_DOMAIN_ERROR
|
|
cout << "BOOST_MATH_THROW_ON_DOMAIN_ERROR" << " is defined to throw on domain error." << endl;
|
|
#else
|
|
cout << "BOOST_MATH_THROW_ON_DOMAIN_ERROR" << " is NOT defined, so NO throw on domain error." << endl;
|
|
#endif
|
|
|
|
// This is a visual sanity check that everything is OK:
|
|
cout << setprecision(17) << showpoint << endl;
|
|
// Show double max_digits10 precision, including trailing zeros.
|
|
|
|
// Test some simple double only examples.
|
|
negative_binomial_distribution<double> my8dist(8., 0.25);
|
|
// 8 successes (r), 0.25 success fraction = 35% or 1 in 4 successes.
|
|
// Note: double values (matching the distribution definition) avoid the need for any casting.
|
|
|
|
BOOST_CHECK_EQUAL(my8dist.successes(), static_cast<double>(8));
|
|
BOOST_CHECK_EQUAL(my8dist.success_fraction(), static_cast<double>(1./4.)); // Exact.
|
|
|
|
double tol = boost::math::tools::epsilon<double>() * 100 * 10;
|
|
// * 100 for %, so tol is 10 epsilon.
|
|
BOOST_CHECK_SMALL(cdf(my8dist, 2.), 4.1580200195313E-4);
|
|
BOOST_CHECK_SMALL(cdf(my8dist, 8.), 0.027129956288264);
|
|
BOOST_CHECK_CLOSE(cdf(my8dist, 16.), 0.233795830683125, tol);
|
|
|
|
BOOST_CHECK_CLOSE(pdf(
|
|
negative_binomial_distribution<double>(2., 0.5),
|
|
1.),
|
|
static_cast<double>(0.25),
|
|
tol);
|
|
|
|
BOOST_CHECK_CLOSE(cdf(complement(
|
|
negative_binomial_distribution<double>(2., 0.5), //
|
|
1.)), // k
|
|
static_cast<double>(0.5), // half
|
|
tol);
|
|
|
|
BOOST_CHECK_CLOSE(cdf(
|
|
negative_binomial_distribution<double>(2., 0.5),
|
|
1.),
|
|
static_cast<double>(0.5),
|
|
tol);
|
|
|
|
BOOST_CHECK_CLOSE(
|
|
quantile(
|
|
negative_binomial_distribution<double>(2., 0.5),
|
|
0.5),
|
|
1., //
|
|
tol);
|
|
|
|
BOOST_CHECK_CLOSE(
|
|
cdf(
|
|
negative_binomial_distribution<double>(8., 0.25),
|
|
16.),
|
|
0.233795830683125, //
|
|
tol);
|
|
|
|
BOOST_CHECK_CLOSE(
|
|
quantile(
|
|
negative_binomial_distribution<double>(8., 0.25),
|
|
0.233795830683125),
|
|
16., //
|
|
tol);
|
|
|
|
// Compare pdf using the simple formulae:
|
|
// exp(lgamma(r + k) - lgamma(r) - lgamma(k+1)) * pow(p, r) * pow((1-p), k)
|
|
// with
|
|
// (p/(r+k) * ibeta_derivative(r, k+1, p) (as used in pdf)
|
|
{
|
|
typedef double RealType;
|
|
RealType r = my8dist.successes();
|
|
RealType p = my8dist.success_fraction();
|
|
for (int i = 0; i <= r * 4; i++) //
|
|
{
|
|
RealType k = static_cast<RealType>(i);
|
|
RealType pmf = exp(lgamma(r + k) - lgamma(r) - lgamma(k+1)) * pow(p, r) * pow((1-p), k);
|
|
BOOST_CHECK_CLOSE(pdf(my8dist, static_cast<RealType>(k)), pmf, tol * 10);
|
|
// 0.0015932321548461931
|
|
// 0.0015932321548461866
|
|
}
|
|
|
|
// Double-check consistency of CDF and PDF by computing the finite sum of pdfs:
|
|
RealType sum = 0;
|
|
for(unsigned i = 0; i <= 20; ++i)
|
|
{
|
|
sum += pdf(my8dist, RealType(i));
|
|
}
|
|
|
|
cout << setprecision(17) << showpoint << sum <<' ' // 0.40025683281803714
|
|
<< cdf(my8dist, static_cast<RealType>(20)) << endl; // 0.40025683281803681
|
|
BOOST_CHECK_CLOSE(sum, cdf(my8dist, static_cast<RealType>(20)), tol);
|
|
}
|
|
|
|
// (Parameter value, arbitrarily zero, only communicates the floating point type).
|
|
// test_spots(0.0F); // Test float.
|
|
test_spots(0.0); // Test double.
|
|
// test_spots(0.0L); // Test long double.
|
|
//#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
|
|
// test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
|
|
//#endif
|
|
|
|
return 0;
|
|
} // int test_main(int, char* [])
|
|
|
|
/*
|
|
|
|
// Test negative binomial cdf using MathCAD computed values.
|
|
// negative_binomial(k, r, p)
|
|
|
|
BOOST_CHECK_CLOSE(negative_binomial_distribution<RealType>(dist(2, 0.25)),
|
|
static_cast<RealType>(0), // k.
|
|
static_cast<RealType>(2), // r
|
|
static_cast<RealType>(0.5)), // mean
|
|
static_cast<RealType>(0.25), // probability 1/4
|
|
tol);
|
|
|
|
BOOST_CHECK_CLOSE(negative_binomial(
|
|
static_cast<RealType>(1), // k
|
|
static_cast<RealType>(4), // r
|
|
static_cast<RealType>(0.5)), // p
|
|
static_cast<RealType>(0.1875), //
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE(negative_binomial(
|
|
static_cast<RealType>(1), // k
|
|
static_cast<RealType>(2), // r
|
|
static_cast<RealType>(0.5)), // pnbinom(1,2,0.5) = 0.5
|
|
static_cast<RealType>(0.5), //
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE(negative_binomial(
|
|
static_cast<RealType>(2), // k
|
|
static_cast<RealType>(2), // n
|
|
static_cast<RealType>(0.5)), // pnbinom(1,2,0.5) = 0.5
|
|
static_cast<RealType>(0.6875), //
|
|
tolerance);
|
|
BOOST_CHECK_CLOSE(negative_binomial(
|
|
static_cast<RealType>(1), // k
|
|
static_cast<RealType>(3), // n
|
|
static_cast<RealType>(0.5)), // pnbinom(1,3,0.5) = 0.5
|
|
static_cast<RealType>(0.3125), // 0.6875
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE(negative_binomial(
|
|
static_cast<RealType>(2), // k
|
|
static_cast<RealType>(2), // n
|
|
static_cast<RealType>(0.1)), //
|
|
static_cast<RealType>(0.0523), //
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE(negative_binomial_c(
|
|
static_cast<RealType>(2), // k
|
|
static_cast<RealType>(2), // n
|
|
static_cast<RealType>(0.1)), //
|
|
static_cast<RealType>(1- 0.0523), //
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE(negative_binomial(
|
|
static_cast<RealType>(26), // k
|
|
static_cast<RealType>(5), // n
|
|
static_cast<RealType>(0.4)), // p
|
|
static_cast<RealType>(0.998968624661119),
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE(negative_binomial(
|
|
static_cast<RealType>(25), // k
|
|
static_cast<RealType>(5), // n
|
|
static_cast<RealType>(0.4)), // p
|
|
static_cast<RealType>(0.998489925933617),
|
|
tolerance);
|
|
|
|
BOOST_CHECK_CLOSE(negative_binomial_c( // complement.
|
|
static_cast<RealType>(26), // k
|
|
static_cast<RealType>(5), // n
|
|
static_cast<RealType>(0.4)), // p
|
|
static_cast<RealType>(1 - 0.998968624661119),
|
|
tolerance * 100); // fails at tolerance * 10
|
|
// 0.0010313753388809799
|
|
*/
|
|
|