123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861 |
- // 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/tools/test.hpp> // for real_concept
- #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
- #define BOOST_TEST_MAIN
- #include <boost/test/unit_test.hpp> // for test_main
- #include <boost/test/tools/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE
- #include "table_type.hpp"
- #include "test_out_of_range.hpp"
- #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 * 100000.0f;
- 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)
- { // 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.
- // Error tests:
- check_out_of_range<negative_binomial_distribution<RealType> >(20, 0.5);
- BOOST_MATH_CHECK_THROW(negative_binomial_distribution<RealType>(0, 0.5), std::domain_error);
- BOOST_MATH_CHECK_THROW(negative_binomial_distribution<RealType>(-2, 0.5), std::domain_error);
- BOOST_MATH_CHECK_THROW(negative_binomial_distribution<RealType>(20, -0.5), std::domain_error);
- BOOST_MATH_CHECK_THROW(negative_binomial_distribution<RealType>(20, 1.5), std::domain_error);
- }
- // 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)),
- static_cast<RealType>(0));
- BOOST_CHECK_EQUAL(
- quantile( // min P < cdf(0) so should be exactly zero.
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(boost::math::tools::min_value<RealType>())),
- static_cast<RealType>(0));
- BOOST_CHECK_CLOSE_FRACTION(
- quantile( // Small P < cdf(0) so should be near zero.
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(boost::math::tools::epsilon<RealType>())), //
- static_cast<RealType>(0),
- tol5eps);
- BOOST_CHECK_CLOSE(
- quantile( // Small 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.0001)),
- static_cast<RealType>(0.95854156929288470),
- tolerance);
- //BOOST_CHECK( // Fails with overflow for real_concept
- //quantile( // Small P near 1 so k failures should be big.
- //negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- //static_cast<RealType>(1 - boost::math::tools::epsilon<RealType>())) <=
- //static_cast<RealType>(189.56999032670058) // 106.462769 for float
- //);
- if(std::numeric_limits<RealType>::has_infinity)
- { // BOOST_CHECK tests for infinity using std::numeric_limits<>::infinity()
- // Note that infinity is not implemented for real_concept, so these tests
- // are only done for types, like built-in float, double.. that have infinity.
- // Note that these assume that BOOST_MATH_OVERFLOW_ERROR_POLICY is NOT throw_on_error.
- // #define BOOST_MATH_THROW_ON_OVERFLOW_POLICY == throw_on_error would throw here.
- // #define BOOST_MAT_DOMAIN_ERROR_POLICY IS defined throw_on_error,
- // so the throw path of error handling is tested below with BOOST_MATH_CHECK_THROW tests.
- BOOST_CHECK(
- quantile( // At P == 1 so k failures should be infinite.
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(1)) ==
- //static_cast<RealType>(boost::math::tools::infinity<RealType>())
- static_cast<RealType>(std::numeric_limits<RealType>::infinity()) );
- BOOST_CHECK_EQUAL(
- quantile( // At 1 == P so should be infinite.
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(1)), //
- std::numeric_limits<RealType>::infinity() );
- BOOST_CHECK_EQUAL(
- quantile(complement( // Q zero 1 so P == 1 < cdf(0) so should be exactly infinity.
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(0))),
- std::numeric_limits<RealType>::infinity() );
- } // test for infinity using std::numeric_limits<>::infinity()
- else
- { // real_concept case, so check it throws rather than returning infinity.
- BOOST_CHECK_EQUAL(
- quantile( // At P == 1 so k failures should be infinite.
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(1)),
- boost::math::tools::max_value<RealType>() );
- BOOST_CHECK_EQUAL(
- quantile(complement( // Q zero 1 so P == 1 < cdf(0) so should be exactly infinity.
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(0))),
- boost::math::tools::max_value<RealType>());
- }
- BOOST_CHECK( // Should work for built-in and real_concept.
- quantile(complement( // Q very near to 1 so P nearly 1 < so should be large > 384.
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(boost::math::tools::min_value<RealType>())))
- >= static_cast<RealType>(384) );
- BOOST_CHECK_EQUAL(
- quantile( // P == 0 < cdf(0) so should be zero.
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(0)),
- static_cast<RealType>(0));
- // Quantile Complement boundary cases:
- BOOST_CHECK_EQUAL(
- quantile(complement( // Q = 1 so P = 0 < cdf(0) so should be exactly zero.
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(1))),
- static_cast<RealType>(0)
- );
- BOOST_CHECK_EQUAL(
- quantile(complement( // Q very near 1 so P == epsilon < cdf(0) so should be exactly zero.
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(1 - boost::math::tools::epsilon<RealType>()))),
- static_cast<RealType>(0)
- );
- // Check that duff arguments throw domain_error:
- BOOST_MATH_CHECK_THROW(
- pdf( // Negative successes!
- negative_binomial_distribution<RealType>(static_cast<RealType>(-1), static_cast<RealType>(0.25)),
- static_cast<RealType>(0)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- pdf( // Negative success_fraction!
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
- static_cast<RealType>(0)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- pdf( // Success_fraction > 1!
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
- static_cast<RealType>(0)),
- std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- pdf( // Negative k argument !
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(-1)),
- std::domain_error
- );
- //BOOST_MATH_CHECK_THROW(
- //pdf( // Unlike binomial there is NO limit on k (failures)
- //negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- //static_cast<RealType>(9)), std::domain_error
- //);
- BOOST_MATH_CHECK_THROW(
- cdf( // Negative k argument !
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
- static_cast<RealType>(-1)),
- std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- cdf( // Negative success_fraction!
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
- static_cast<RealType>(0)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- cdf( // Success_fraction > 1!
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
- static_cast<RealType>(0)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- quantile( // Negative success_fraction!
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(-0.25)),
- static_cast<RealType>(0)), std::domain_error
- );
- BOOST_MATH_CHECK_THROW(
- quantile( // Success_fraction > 1!
- negative_binomial_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(1.25)),
- static_cast<RealType>(0)), std::domain_error
- );
- // End of check throwing 'duff' out-of-domain values.
- #define T RealType
- #include "negative_binomial_quantile.ipp"
- for(unsigned i = 0; i < negative_binomial_quantile_data.size(); ++i)
- {
- using namespace boost::math::policies;
- typedef policy<discrete_quantile<boost::math::policies::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>() * 700;
- if(!boost::is_floating_point<RealType>::value)
- tol *= 10; // no lanczos approximation implies less accuracy
- //
- // Check full real value first:
- //
- negative_binomial_distribution<RealType, P1> p1(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]);
- RealType x = quantile(p1, negative_binomial_quantile_data[i][2]);
- BOOST_CHECK_CLOSE_FRACTION(x, negative_binomial_quantile_data[i][3], tol);
- x = quantile(complement(p1, negative_binomial_quantile_data[i][2]));
- BOOST_CHECK_CLOSE_FRACTION(x, negative_binomial_quantile_data[i][4], tol);
- //
- // Now with round down to integer:
- //
- negative_binomial_distribution<RealType, P2> p2(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]);
- x = quantile(p2, negative_binomial_quantile_data[i][2]);
- BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][3]));
- x = quantile(complement(p2, negative_binomial_quantile_data[i][2]));
- BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][4]));
- //
- // Now with round up to integer:
- //
- negative_binomial_distribution<RealType, P3> p3(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]);
- x = quantile(p3, negative_binomial_quantile_data[i][2]);
- BOOST_CHECK_EQUAL(x, ceil(negative_binomial_quantile_data[i][3]));
- x = quantile(complement(p3, negative_binomial_quantile_data[i][2]));
- BOOST_CHECK_EQUAL(x, ceil(negative_binomial_quantile_data[i][4]));
- //
- // Now with round to integer "outside":
- //
- negative_binomial_distribution<RealType, P4> p4(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]);
- x = quantile(p4, negative_binomial_quantile_data[i][2]);
- BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? floor(negative_binomial_quantile_data[i][3]) : ceil(negative_binomial_quantile_data[i][3]));
- x = quantile(complement(p4, negative_binomial_quantile_data[i][2]));
- BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? ceil(negative_binomial_quantile_data[i][4]) : floor(negative_binomial_quantile_data[i][4]));
- //
- // Now with round to integer "inside":
- //
- negative_binomial_distribution<RealType, P5> p5(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]);
- x = quantile(p5, negative_binomial_quantile_data[i][2]);
- BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? ceil(negative_binomial_quantile_data[i][3]) : floor(negative_binomial_quantile_data[i][3]));
- x = quantile(complement(p5, negative_binomial_quantile_data[i][2]));
- BOOST_CHECK_EQUAL(x, negative_binomial_quantile_data[i][2] < 0.5f ? floor(negative_binomial_quantile_data[i][4]) : ceil(negative_binomial_quantile_data[i][4]));
- //
- // Now with round to nearest integer:
- //
- negative_binomial_distribution<RealType, P6> p6(negative_binomial_quantile_data[i][0], negative_binomial_quantile_data[i][1]);
- x = quantile(p6, negative_binomial_quantile_data[i][2]);
- BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][3] + 0.5f));
- x = quantile(complement(p6, negative_binomial_quantile_data[i][2]));
- BOOST_CHECK_EQUAL(x, floor(negative_binomial_quantile_data[i][4] + 0.5f));
- }
- return;
- } // template <class RealType> void test_spots(RealType) // Any floating-point type RealType.
- BOOST_AUTO_TEST_CASE( test_main )
- {
- // 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.
- // 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.
- // Check accessor functions return exact values for double at least.
- BOOST_CHECK_EQUAL(my8dist.successes(), static_cast<double>(8));
- BOOST_CHECK_EQUAL(my8dist.success_fraction(), static_cast<double>(1./4.));
- // (Parameter value, arbitrarily zero, only communicates the floating point type).
- #ifdef TEST_FLOAT
- 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
- #ifdef TEST_LDOUBLE
- test_spots(0.0L); // Test long double.
- #endif
- #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
- #ifdef TEST_REAL_CONCEPT
- test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
- #endif
- #endif
- #else
- std::cout << "<note>The long double tests have been disabled on this platform "
- "either because the long double overloads of the usual math functions are "
- "not available at all, or because they are too inaccurate for these tests "
- "to pass.</note>" << std::endl;
- #endif
-
- } // BOOST_AUTO_TEST_CASE( test_main )
- /*
- Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_negative_binomial.exe"
- Running 1 test case...
- Tolerance = 0.0119209%.
- Tolerance 5 eps = 5.96046e-007%.
- Tolerance = 2.22045e-011%.
- Tolerance 5 eps = 1.11022e-015%.
- Tolerance = 2.22045e-011%.
- Tolerance 5 eps = 1.11022e-015%.
- Tolerance = 2.22045e-011%.
- Tolerance 5 eps = 1.11022e-015%.
- *** No errors detected
- */
|