diff options
Diffstat (limited to 'src/boost/libs/math/test/test_geometric.cpp')
| -rw-r--r-- | src/boost/libs/math/test/test_geometric.cpp | 802 |
1 files changed, 802 insertions, 0 deletions
diff --git a/src/boost/libs/math/test/test_geometric.cpp b/src/boost/libs/math/test/test_geometric.cpp new file mode 100644 index 00000000..33fbb8d4 --- /dev/null +++ b/src/boost/libs/math/test/test_geometric.cpp @@ -0,0 +1,802 @@ +// test_geometric.cpp + +// Copyright Paul A. Bristow 2010. +// Copyright John Maddock 2010. + +// 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 Geometric 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> +#include <boost/math/concepts/real_concept.hpp> // for real_concept +using ::boost::math::concepts::real_concept; + +#include <boost/math/distributions/geometric.hpp> // for geometric_distribution +using boost::math::geometric_distribution; +using boost::math::geometric; // using typedef for geometric_distribution<double> + +#include <boost/math/distributions/negative_binomial.hpp> // for some comparisons. + +#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_FRACTION +#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 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::geometric_distribution<RealType> g(p); + BOOST_CHECK_EQUAL(p, g.success_fraction()); + BOOST_CHECK_CLOSE_FRACTION(cdf(g, 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_FRACTION( + cdf(complement(g, k)), Q, tol); + if(k != 0) + { + BOOST_CHECK_CLOSE_FRACTION( + quantile(g, 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(g, P) < boost::math::tools::epsilon<RealType>() * 10); + } + } + if(k != 0) + { + BOOST_CHECK_CLOSE_FRACTION( + quantile(complement(g, 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(g, Q)) < boost::math::tools::epsilon<RealType>() * 10); + } + } + } // if((P < 0.99) && (Q < 0.99)) + + // Parameter estimation test: estimate success ratio: + BOOST_CHECK_CLOSE_FRACTION( + geometric_distribution<RealType>::find_lower_bound_on_p( + 1+k, P), + p, 0.02); // Wide tolerance needed for some tests. + // 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_FRACTION( + geometric_distribution<RealType>::find_upper_bound_on_p( + 1+k+1, Q), + p, 0.02); + + 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( + geometric_distribution<RealType>::find_lower_bound_on_p( + 1+k, Q) + <= + geometric_distribution<RealType>::find_upper_bound_on_p( + 1+k, Q) + ); + BOOST_CHECK( + geometric_distribution<RealType>::find_lower_bound_on_p( + 1+k, Q) + <= + 1 / (1+k) + ); + BOOST_CHECK( + 1 / (1+k) + <= + geometric_distribution<RealType>::find_upper_bound_on_p( + 1+k, Q) + ); + } + else + { + // As above but when P is small. + BOOST_CHECK( + geometric_distribution<RealType>::find_lower_bound_on_p( + 1+k, P) + <= + geometric_distribution<RealType>::find_upper_bound_on_p( + 1+k, P) + ); + BOOST_CHECK( + geometric_distribution<RealType>::find_lower_bound_on_p( + 1+k, P) + <= + 1 / (1+k) + ); + BOOST_CHECK( + 1 / (1+k) + <= + geometric_distribution<RealType>::find_upper_bound_on_p( + 1+k, P) + ); + } + + // Estimate sample size: + BOOST_CHECK_CLOSE_FRACTION( + geometric_distribution<RealType>::find_minimum_number_of_trials( + k, p, P), + 1+k, 0.02); // Can differ 50 to 51 for small p + BOOST_CHECK_CLOSE_FRACTION( + geometric_distribution<RealType>::find_maximum_number_of_trials( + k, p, Q), + 1+k, 0.02); + +} // test_spot + +template <class RealType> // Any floating-point type RealType. +void test_spots(RealType) +{ + // Basic sanity checks. + // Most test data is to double precision (17 decimal digits) only, + + cout << "Floating point Type is " << typeid(RealType).name() << endl; + + // so set tolerance to 1000 eps expressed as a fraction, + // or 1000 eps of type double expressed as a fraction, + // whichever is the larger. + + RealType tolerance = (std::max) + (boost::math::tools::epsilon<RealType>(), + static_cast<RealType>(std::numeric_limits<double>::epsilon())); + tolerance *= 10; // 10 eps + + cout << "Tolerance = " << tolerance << "." << endl; + + RealType tol1eps = boost::math::tools::epsilon<RealType>(); // Very tight, suit exact values. + //RealType tol2eps = boost::math::tools::epsilon<RealType>() * 2; // Tight, values. + RealType tol5eps = boost::math::tools::epsilon<RealType>() * 5; // Wider 5 epsilon. + cout << "Tolerance 5 eps = " << tol5eps << "." << endl; + + + // Sources of spot test values are mainly R. + + using boost::math::geometric_distribution; + using boost::math::geometric; + using boost::math::cdf; + using boost::math::pdf; + using boost::math::quantile; + using boost::math::complement; + + BOOST_MATH_STD_USING // for std math functions + + // Test geometric using cdf spot values R + // These test quantiles and complements as well. + + test_spot( // + static_cast<RealType>(2), // Number of failures, k + static_cast<RealType>(0.5), // Probability of success as fraction, p + static_cast<RealType>(0.875L), // Probability of result (CDF), P + static_cast<RealType>(0.125L), // complement CCDF Q = 1 - P + tolerance); + + test_spot( // + static_cast<RealType>(0), // Number of failures, k + static_cast<RealType>(0.25), // Probability of success as fraction, p + static_cast<RealType>(0.25), // Probability of result (CDF), P + static_cast<RealType>(0.75), // Q = 1 - P + tolerance); + + test_spot( + // R formatC(pgeom(10,0.25), digits=17) [1] "0.95776486396789551" + // formatC(pgeom(10,0.25, FALSE), digits=17) [1] "0.042235136032104499" + + static_cast<RealType>(10), // Number of failures, k + static_cast<RealType>(0.25), // Probability of success, p + static_cast<RealType>(0.95776486396789551L), // Probability of result (CDF), P + static_cast<RealType>(0.042235136032104499L), // Q = 1 - P + tolerance); + + test_spot( // + // > R formatC(pgeom(50,0.25, TRUE), digits=17) [1] "0.99999957525875771" + // > R formatC(pgeom(50,0.25, FALSE), digits=17) [1] "4.2474124232020353e-07" + static_cast<RealType>(50), // Number of failures, k + static_cast<RealType>(0.25), // Probability of success, p + static_cast<RealType>(0.99999957525875771), // Probability of result (CDF), P + static_cast<RealType>(4.2474124232020353e-07), // Q = 1 - P + tolerance); + /* + // This causes failures in find_upper_bound_on_p p is small branch. + test_spot( // formatC(pgeom(50,0.01, TRUE), digits=17)[1] "0.40104399353383874" + // > formatC(pgeom(50,0.01, FALSE), digits=17) [1] "0.59895600646616121" + static_cast<RealType>(50), // Number of failures, k + static_cast<RealType>(0.01), // Probability of success, p + static_cast<RealType>(0.40104399353383874), // Probability of result (CDF), P + static_cast<RealType>(0.59895600646616121), // Q = 1 - P + tolerance); + */ + + test_spot( // > formatC(pgeom(50,0.99, TRUE), digits=17) [1] " 1" + // formatC(pgeom(50,0.99, FALSE), digits=17) [1] "1.0000000000000364e-102" + static_cast<RealType>(50), // Number of failures, k + static_cast<RealType>(0.99), // Probability of success, p + static_cast<RealType>(1), // Probability of result (CDF), P + static_cast<RealType>(1.0000000000000364e-102), // Q = 1 - P + tolerance); + + test_spot( // > formatC(pgeom(1,0.99, TRUE), digits=17) [1] "0.99990000000000001" + // > formatC(pgeom(1,0.99, FALSE), digits=17) [1] "0.00010000000000000009" + static_cast<RealType>(1), // Number of failures, k + static_cast<RealType>(0.99), // Probability of success, p + static_cast<RealType>(0.9999), // Probability of result (CDF), P + static_cast<RealType>(0.0001), // Q = 1 - P + tolerance); + +if(std::numeric_limits<RealType>::is_specialized) +{ // An extreme value test that is more accurate than using negative binomial. + // Since geometric only uses exp and log functions. + test_spot( // > formatC(pgeom(10000, 0.001, TRUE), digits=17) [1] "0.99995487182736897" +// > formatC(pgeom(10000,0.001, FALSE), digits=17) [1] "4.5128172631071587e-05" + static_cast<RealType>(10000L), // Number of failures, k + static_cast<RealType>(0.001L), // Probability of success, p + static_cast<RealType>(0.99995487182736897L), // Probability of result (CDF), P + static_cast<RealType>(4.5128172631071587e-05L), // Q = 1 - P + tolerance); // + } // numeric_limit is specialized + // End of single spot tests using RealType + + // Tests on PDF: + + BOOST_CHECK_CLOSE_FRACTION( //> formatC(dgeom(0,0.5), digits=17)[1] " 0.5" + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)), + static_cast<RealType>(0.0) ), // Number of failures, k is very small but not integral, + static_cast<RealType>(0.5), // nearly success probability. + tolerance); + + BOOST_CHECK_CLOSE_FRACTION( //> formatC(dgeom(0,0.5), digits=17)[1] " 0.5" + // R treates geom as a discrete distribution. + // > formatC(dgeom(1.999999,0.5, FALSE), digits=17) [1] " 0" + // Warning message: + // In dgeom(1.999999, 0.5, FALSE) : non-integer x = 1.999999 + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)), + static_cast<RealType>(0.0001L) ), // Number of failures, k is very small but not integral, + static_cast<RealType>(0.4999653438420768L), // nearly success probability. + tolerance); + + BOOST_CHECK_CLOSE_FRACTION( // > formatC(pgeom(0.0001,0.5, TRUE), digits=17)[1] " 0.5" + // > formatC(pgeom(0.0001,0.5, FALSE), digits=17) [1] " 0.5" + // R treates geom as a discrete distribution. + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)), + static_cast<RealType>(0.0001L) ), // Number of failures, k is very small but not integral, + static_cast<RealType>(0.4999653438420768L), // nearly success probability. + tolerance); + + BOOST_CHECK_CLOSE_FRACTION( // formatC(dgeom(1,0.01), digits=17)[1] "0.0099000000000000008" + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.01L)), + static_cast<RealType>(1) ), // Number of failures, k + static_cast<RealType>(0.0099000000000000008), // + tolerance); + + BOOST_CHECK_CLOSE_FRACTION( //> formatC(dgeom(1,0.99), digits=17)[1] "0.0099000000000000043" + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.99L)), + static_cast<RealType>(1) ), // Number of failures, k + static_cast<RealType>(0.00990000000000000043L), // + tolerance); + + BOOST_CHECK_CLOSE_FRACTION( //> > formatC(dgeom(0,0.99), digits=17)[1] "0.98999999999999999" + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.99L)), + static_cast<RealType>(0) ), // Number of failures, k + static_cast<RealType>(0.98999999999999999L), // + tolerance); + + // p near unity. + BOOST_CHECK_CLOSE_FRACTION( // > formatC(dgeom(100,0.99), digits=17)[1] "9.9000000000003448e-201" + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.99L)), + static_cast<RealType>(100) ), // Number of failures, k + static_cast<RealType>(9.9000000000003448e-201L), // + 100 * tolerance); // Note difference + + // p nearer unity. + BOOST_CHECK_CLOSE_FRACTION( // + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999)), + static_cast<RealType>(10) ), // Number of failures, k + // static_cast<double>(9.9989999999889024e-41), // Boost.Math + // static_cast<float>(1.00156406e-040) + static_cast<RealType>(9.999e-41), // exact from 100 digit calculator. + 2e3 * tolerance); // Note bigger tolerance needed. + + // Moshier Cephes 100 digits calculator says 9.999e-41 + //0.9999*pow(1-0.9999,10) + // 9.9990000000000000000000000000000000000000000000000000000000000000000000E-41 + // 9.998999999988988e-041 + // > formatC(dgeom(10, 0.9999), digits=17) [1] "9.9989999999889024e-41" + // p * pow(q, k) 9.9989999999889880e-041 + // exp(p * k * log1p(-p)) 9.9989999999889024e-041 + + + + // 0.9999999999 * pow(1-0.9999999999,10)= 9.9999999990E-101 + // > formatC(dgeom(10,0.9999999999), digits=17) [1] "1.0000008273040127e-100" + BOOST_CHECK_CLOSE_FRACTION( // + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999999999L)), + static_cast<RealType>(10) ), // + static_cast<RealType>(9.9999999990E-101L), // 1.0000008273040179e-100 + 1e9 * tolerance); // Note big tolerance needed. + // 1.0000008273040179e-100 Boost.Math + // 1.0000008273040127e-100 R + // 0.9999999990000004e-100 100 digit calculator 'exact' + + BOOST_CHECK_CLOSE_FRACTION( // + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.00000000001L)), + static_cast<RealType>(10) ), // + static_cast<RealType>(9.999999999e-12L), // get 9.9999999989999994e-012 + 1 * tolerance); // Note small tolerance needed. + + + BOOST_CHECK_CLOSE_FRACTION( // + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.00000000001L)), + static_cast<RealType>(1000) ), // + static_cast<RealType>(9.9999999e-12L), // get 9.9999998999999913e-012 + tolerance); // Note small tolerance needed. + + + /////////////////////////////////////////////////// + BOOST_CHECK_CLOSE_FRACTION( // + // > formatC(dgeom(0.0001,0.5, FALSE), digits=17) [1] " 0.5" + // R treates geom as a discrete distribution. + // But Boost.Math is continuous, so if you want R behaviour, + // make number of failures, k into an integer with the floor function. + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)), + static_cast<RealType>(floor(0.0001L)) ), // Number of failures, k is very small but MADE integral, + static_cast<RealType>(0.5), // nearly success probability. + tolerance); + + // R switches over at about 1e7 from k = 0, returning 0.5, to k = 1, returning 0.25. + // Boost.Math does not do this, even for 0.9999999999999999 + // > formatC(pgeom(0.999999,0.5, FALSE), digits=17) [1] " 0.5" + // > formatC(pgeom(0.9999999,0.5, FALSE), digits=17) [1] " 0.25" + + BOOST_CHECK_CLOSE_FRACTION( // > formatC(pgeom(0.0001,0.5, TRUE), digits=17)[1] " 0.5" + // > formatC(pgeom(0.0001,0.5, FALSE), digits=17) [1] " 0.5" + // R treates geom as a discrete distribution. + // But Boost.Math is continuous, so if you want R behaviour, + // make number of failures, k into an integer with the floor function. + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)), + static_cast<RealType>(floor(0.9999999999999999L)) ), // Number of failures, k is very small but MADE integral, + static_cast<RealType>(0.5), // nearly success probability. + tolerance); + + BOOST_CHECK_CLOSE_FRACTION( // > formatC(pgeom(0.0001,0.5, TRUE), digits=17)[1] " 0.5" + // > formatC(pgeom(0.0001,0.5, FALSE), digits=17) [1] " 0.5" + // R treates geom as a discrete distribution. + // But Boost.Math is continuous, so if you want R behaviour, + // make number of failures, k into an integer with the floor function. + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)), + static_cast<RealType>(floor(1. - tolerance)) ), + // Number of failures, k is very small but MADE integral, + // Need to use tolerance here, + // as epsilon is ill-defined for Real concept: + // numeric_limits<RealType>::epsilon() 0 + static_cast<RealType>(0.5), // nearly success probability. + tolerance * 10); + + BOOST_CHECK_CLOSE_FRACTION( + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.0001L)), + static_cast<RealType>(2)), // k = 2. + static_cast<RealType>(9.99800010e-5L), // 'exact ' + tolerance); + + //> formatC(dgeom(2, 0.9999), digits=17) [1] "9.9989999999977806e-09" + BOOST_CHECK_CLOSE_FRACTION( + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999L)), + static_cast<RealType>(2)), // k = 0 + static_cast<RealType>(9.999e-9L), // 'exact' + 1000*tolerance); + + BOOST_CHECK_CLOSE_FRACTION( + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999L)), + static_cast<RealType>(3)), // k = 3 + static_cast<RealType>(9.999e-13L), // get + 1000*tolerance); + + BOOST_CHECK_CLOSE_FRACTION( + pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999L)), + static_cast<RealType>(5)), // k = 5 + static_cast<RealType>(9.999e-21L), // 9.9989999999944947e-021 + 1000*tolerance); + + + BOOST_CHECK_CLOSE_FRACTION( + pdf(geometric_distribution<RealType>( static_cast<RealType>(0.0001L)), + static_cast<RealType>(3)), // k = 0. + static_cast<RealType>(9.99700029999e-5L), // + tolerance); + // Tests on cdf: + // MathCAD pgeom k, r, p) == failures, successes, probability. + + BOOST_CHECK_CLOSE_FRACTION(cdf( + geometric_distribution<RealType>(static_cast<RealType>(0.5)), // prob 0.5 + static_cast<RealType>(0) ), // k = 0 + static_cast<RealType>(0.5), // probability =p + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(cdf(complement( + geometric_distribution<RealType>(static_cast<RealType>(0.5)), // + static_cast<RealType>(0) )), // k = 0 + static_cast<RealType>(0.5), // probability = + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(cdf( + geometric_distribution<RealType>(static_cast<RealType>(0.25)), // prob 0.5 + static_cast<RealType>(1) ), // k = 0 + static_cast<RealType>(0.4375L), // probability =p + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(cdf(complement( + geometric_distribution<RealType>(static_cast<RealType>(0.25)), // + static_cast<RealType>(1) )), // k = 0 + static_cast<RealType>(1-0.4375L), // probability = + tolerance); + + BOOST_CHECK_CLOSE_FRACTION(cdf(complement( + geometric_distribution<RealType>(static_cast<RealType>(0.5)), // + static_cast<RealType>(1) )), // k = 0 + static_cast<RealType>(0.25), // probability = exact 0.25 + tolerance); + + BOOST_CHECK_CLOSE_FRACTION( // + cdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)), + static_cast<RealType>(4)), // k =4. + static_cast<RealType>(0.96875L), // exact + tolerance); + + + // Tests of other functions, mean and other moments ... + + geometric_distribution<RealType> dist(static_cast<RealType>(0.25)); + // mean: + BOOST_CHECK_CLOSE_FRACTION( + mean(dist), static_cast<RealType>((1 - 0.25) /0.25), tol5eps); + BOOST_CHECK_CLOSE_FRACTION( + mode(dist), static_cast<RealType>(0), tol1eps); + // variance: + BOOST_CHECK_CLOSE_FRACTION( + variance(dist), static_cast<RealType>((1 - 0.25) / (0.25 * 0.25)), tol5eps); + + // std deviation: + // sqrt(0.75/0.125) + + BOOST_CHECK_CLOSE_FRACTION( + standard_deviation(dist), // + static_cast<RealType>(sqrt((1.0L - 0.25L) / (0.25L * 0.25L))), // using 100 digit calc + tol5eps); + + BOOST_CHECK_CLOSE_FRACTION( + skewness(dist), // + static_cast<RealType>((2-0.25L) /sqrt(0.75L)), + // using calculator + tol5eps); + BOOST_CHECK_CLOSE_FRACTION( + kurtosis_excess(dist), // + static_cast<RealType>(6 + 0.0625L/0.75L), // + tol5eps); + // 6.083333333333333 6.166666666666667 + BOOST_CHECK_CLOSE_FRACTION( + kurtosis(dist), // true + static_cast<RealType>(9 + 0.0625L/0.75L), // + tol5eps); + // hazard: + RealType x = static_cast<RealType>(0.125); + BOOST_CHECK_CLOSE_FRACTION( + hazard(dist, x) + , pdf(dist, x) / cdf(complement(dist, x)), tol5eps); + // cumulative hazard: + BOOST_CHECK_CLOSE_FRACTION( + chf(dist, x), -log(cdf(complement(dist, x))), tol5eps); + // coefficient_of_variation: + BOOST_CHECK_CLOSE_FRACTION( + coefficient_of_variation(dist) + , standard_deviation(dist) / mean(dist), tol5eps); + + // Special cases for PDF: + BOOST_CHECK_EQUAL( + pdf( + geometric_distribution<RealType>(static_cast<RealType>(0)), // + static_cast<RealType>(0)), + static_cast<RealType>(0) ); + + BOOST_CHECK_EQUAL( + pdf( + geometric_distribution<RealType>(static_cast<RealType>(0)), + static_cast<RealType>(0.0001)), + static_cast<RealType>(0) ); + + BOOST_CHECK_EQUAL( + pdf( + geometric_distribution<RealType>(static_cast<RealType>(1)), + static_cast<RealType>(0.001)), + static_cast<RealType>(0) ); + + BOOST_CHECK_EQUAL( + pdf( + geometric_distribution<RealType>(static_cast<RealType>(1)), + static_cast<RealType>(8)), + static_cast<RealType>(0) ); + + BOOST_CHECK_SMALL( + pdf( + geometric_distribution<RealType>(static_cast<RealType>(0.25)), + static_cast<RealType>(0))- + static_cast<RealType>(0.25), + 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. + geometric_distribution<RealType>(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. + geometric_distribution<RealType>(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. + geometric_distribution<RealType>(static_cast<RealType>(0.25)), + static_cast<RealType>(boost::math::tools::epsilon<RealType>())), // + static_cast<RealType>(0), + tol5eps); + + BOOST_CHECK_CLOSE_FRACTION( + quantile( // Small P < cdf(0) so should be exactly zero. + geometric_distribution<RealType>(static_cast<RealType>(0.25)), + static_cast<RealType>(0.0001)), + static_cast<RealType>(0), + tolerance); + + //BOOST_CHECK( // Fails with overflow for real_concept + //quantile( // Small P near 1 so k failures should be big. + //geometric_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. + geometric_distribution<RealType>(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. + geometric_distribution<RealType>( 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. + geometric_distribution<RealType>(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. + geometric_distribution<RealType>(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. + geometric_distribution<RealType>(static_cast<RealType>(0.25)), + static_cast<RealType>(0))), + boost::math::tools::max_value<RealType>()); + } // has infinity + + BOOST_CHECK( // Should work for built-in and real_concept. + quantile(complement( // Q near to 1 so P nearly 1, so should be large > 300. + geometric_distribution<RealType>(static_cast<RealType>(0.25)), + static_cast<RealType>(boost::math::tools::min_value<RealType>()))) + >= static_cast<RealType>(300) ); + + BOOST_CHECK_EQUAL( + quantile( // P == 0 < cdf(0) so should be zero. + geometric_distribution<RealType>(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. + geometric_distribution<RealType>( 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. + geometric_distribution<RealType>(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 success_fraction! + geometric_distribution<RealType>(static_cast<RealType>(-0.25)), + static_cast<RealType>(0)), std::domain_error); + BOOST_MATH_CHECK_THROW( + pdf( // Success_fraction > 1! + geometric_distribution<RealType>(static_cast<RealType>(1.25)), + static_cast<RealType>(0)), + std::domain_error); + BOOST_MATH_CHECK_THROW( + pdf( // Negative k argument ! + geometric_distribution<RealType>(static_cast<RealType>(0.25)), + static_cast<RealType>(-1)), + std::domain_error); + //BOOST_MATH_CHECK_THROW( + //pdf( // check limit on k (failures) + //geometric_distribution<RealType>(static_cast<RealType>(0.25)), + //std::numeric_limits<RealType>infinity()), + //std::domain_error); + BOOST_MATH_CHECK_THROW( + cdf( // Negative k argument ! + geometric_distribution<RealType>(static_cast<RealType>(0.25)), + static_cast<RealType>(-1)), + std::domain_error); + BOOST_MATH_CHECK_THROW( + cdf( // Negative success_fraction! + geometric_distribution<RealType>(static_cast<RealType>(-0.25)), + static_cast<RealType>(0)), std::domain_error); + BOOST_MATH_CHECK_THROW( + cdf( // Success_fraction > 1! + geometric_distribution<RealType>(static_cast<RealType>(1.25)), + static_cast<RealType>(0)), std::domain_error); + BOOST_MATH_CHECK_THROW( + quantile( // Negative success_fraction! + geometric_distribution<RealType>(static_cast<RealType>(-0.25)), + static_cast<RealType>(0)), std::domain_error); + BOOST_MATH_CHECK_THROW( + quantile( // Success_fraction > 1! + geometric_distribution<RealType>(static_cast<RealType>(1.25)), + static_cast<RealType>(0)), std::domain_error); + check_out_of_range<geometric_distribution<RealType> >(0.5); + // End of check throwing 'duff' out-of-domain values. + + { // Compare geometric and negative binomial functions. + using boost::math::negative_binomial_distribution; + using boost::math::geometric_distribution; + + RealType k = static_cast<RealType>(2.L); + RealType alpha = static_cast<RealType>(0.05L); + RealType p = static_cast<RealType>(0.5L); + + BOOST_CHECK_CLOSE_FRACTION( // Successes parameter in negative binomial is 1 for geometric. + geometric_distribution<RealType>::find_lower_bound_on_p(k, alpha), + negative_binomial_distribution<RealType>::find_lower_bound_on_p(k, static_cast<RealType>(1), alpha), + tolerance); + BOOST_CHECK_CLOSE_FRACTION( // Successes parameter in negative binomial is 1 for geometric. + geometric_distribution<RealType>::find_upper_bound_on_p(k, alpha), + negative_binomial_distribution<RealType>::find_upper_bound_on_p(k, static_cast<RealType>(1), alpha), + tolerance); + BOOST_CHECK_CLOSE_FRACTION( // Should be identical - successes parameter is not used. + geometric_distribution<RealType>::find_maximum_number_of_trials(k, p, alpha), + negative_binomial_distribution<RealType>::find_maximum_number_of_trials(k, p, alpha), + tolerance); + } + //geometric::find_upper_bound_on_p(k, alpha); + return; +} // template <class RealType> void test_spots(RealType) // Any floating-point type RealType. + +BOOST_AUTO_TEST_CASE( test_main ) +{ + // Check that can generate geometric distribution using the two convenience methods: + using namespace boost::math; + geometric g05d(0.5); // Using typedef - default type is double. + geometric_distribution<> g05dd(0.5); // Using default RealType double. + + // Basic sanity-check spot values. + + // Test some simple double only examples. + geometric_distribution<double> mydist(0.25); + // success fraction == 0.25 == 25% 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(mydist.success_fraction(), static_cast<double>(1./4.)); + + //cout << numeric_limits<RealType>::epsilon() << endl; + + // (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 + #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) +#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 ) + +/* + + + +*/ |