123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444 |
- // (C) Copyright John Maddock 2007.
- // 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)
- #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
- #include <boost/math/concepts/real_concept.hpp>
- #define BOOST_TEST_MAIN
- #include <boost/test/unit_test.hpp>
- #include <boost/test/tools/floating_point_comparison.hpp>
- #include <boost/math/distributions/non_central_chi_squared.hpp>
- #include <boost/type_traits/is_floating_point.hpp>
- #include <boost/array.hpp>
- #include "functor.hpp"
- #include "handle_test_result.hpp"
- #include "table_type.hpp"
- #include <iostream>
- #include <iomanip>
- #define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \
- {\
- unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
- BOOST_CHECK_CLOSE(a, b, prec); \
- if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
- {\
- std::cerr << "Failure was at row " << i << std::endl;\
- std::cerr << std::setprecision(35); \
- std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
- std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
- }\
- }
- #define BOOST_CHECK_EX(a, i) \
- {\
- unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
- BOOST_CHECK(a); \
- if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
- {\
- std::cerr << "Failure was at row " << i << std::endl;\
- std::cerr << std::setprecision(35); \
- std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
- std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
- }\
- }
- template <class RealType>
- RealType naive_pdf(RealType v, RealType lam, RealType x)
- {
- // Formula direct from
- // http://mathworld.wolfram.com/NoncentralChi-SquaredDistribution.html
- // with no simplification:
- RealType sum, term, prefix(1);
- RealType eps = boost::math::tools::epsilon<RealType>();
- term = sum = pdf(boost::math::chi_squared_distribution<RealType>(v), x);
- for(int i = 1;; ++i)
- {
- prefix *= lam / (2 * i);
- term = prefix * pdf(boost::math::chi_squared_distribution<RealType>(v + 2 * i), x);
- sum += term;
- if(term / sum < eps)
- break;
- }
- return sum * exp(-lam / 2);
- }
- template <class RealType>
- void test_spot(
- RealType df, // Degrees of freedom
- RealType ncp, // non-centrality param
- RealType cs, // Chi Square statistic
- RealType P, // CDF
- RealType Q, // Complement of CDF
- RealType tol) // Test tolerance
- {
- boost::math::non_central_chi_squared_distribution<RealType> dist(df, ncp);
- BOOST_CHECK_CLOSE(
- cdf(dist, cs), P, tol);
- #ifndef BOOST_NO_EXCEPTIONS
- try{
- BOOST_CHECK_CLOSE(
- pdf(dist, cs), naive_pdf(dist.degrees_of_freedom(), ncp, cs), tol * 150);
- }
- catch(const std::overflow_error&)
- {
- }
- #endif
- 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 Q is reasonably free of error:
- //
- BOOST_CHECK_CLOSE(
- cdf(complement(dist, cs)), Q, tol);
- BOOST_CHECK_CLOSE(
- quantile(dist, P), cs, tol * 10);
- BOOST_CHECK_CLOSE(
- quantile(complement(dist, Q)), cs, tol * 10);
- BOOST_CHECK_CLOSE(
- dist.find_degrees_of_freedom(ncp, cs, P), df, tol * 10);
- BOOST_CHECK_CLOSE(
- dist.find_degrees_of_freedom(boost::math::complement(ncp, cs, Q)), df, tol * 10);
- BOOST_CHECK_CLOSE(
- dist.find_non_centrality(df, cs, P), ncp, tol * 10);
- BOOST_CHECK_CLOSE(
- dist.find_non_centrality(boost::math::complement(df, cs, Q)), ncp, tol * 10);
- }
- }
- template <class RealType> // Any floating-point type RealType.
- void test_spots(RealType)
- {
- #ifndef ERROR_REPORTING_MODE
- RealType tolerance = (std::max)(
- boost::math::tools::epsilon<RealType>(),
- (RealType)boost::math::tools::epsilon<double>() * 5) * 150;
- //
- // At float precision we need to up the tolerance, since
- // the input values are rounded off to inexact quantities
- // the results get thrown off by a noticeable amount.
- //
- if(boost::math::tools::digits<RealType>() < 50)
- tolerance *= 50;
- if(boost::is_floating_point<RealType>::value != 1)
- tolerance *= 20; // real_concept special functions are less accurate
- std::cout << "Tolerance = " << tolerance << "%." << std::endl;
- using boost::math::chi_squared_distribution;
- using ::boost::math::chi_squared;
- using ::boost::math::cdf;
- using ::boost::math::pdf;
- //
- // Test against the data from Table 6 of:
- //
- // "Self-Validating Computations of Probabilities for Selected
- // Central and Noncentral Univariate Probability Functions."
- // Morgan C. Wang; William J. Kennedy
- // Journal of the American Statistical Association,
- // Vol. 89, No. 427. (Sep., 1994), pp. 878-887.
- //
- test_spot(
- static_cast<RealType>(1), // degrees of freedom
- static_cast<RealType>(6), // non centrality
- static_cast<RealType>(0.00393), // Chi Squared statistic
- static_cast<RealType>(0.2498463724258039e-2), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.2498463724258039e-2), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(5), // degrees of freedom
- static_cast<RealType>(1), // non centrality
- static_cast<RealType>(9.23636), // Chi Squared statistic
- static_cast<RealType>(0.8272918751175548), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.8272918751175548), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(11), // degrees of freedom
- static_cast<RealType>(21), // non centrality
- static_cast<RealType>(24.72497), // Chi Squared statistic
- static_cast<RealType>(0.2539481822183126), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.2539481822183126), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(31), // degrees of freedom
- static_cast<RealType>(6), // non centrality
- static_cast<RealType>(44.98534), // Chi Squared statistic
- static_cast<RealType>(0.8125198785064969), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.8125198785064969), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(51), // degrees of freedom
- static_cast<RealType>(1), // non centrality
- static_cast<RealType>(38.56038), // Chi Squared statistic
- static_cast<RealType>(0.8519497361859118e-1), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.8519497361859118e-1), // Q = 1 - P
- tolerance * 2);
- test_spot(
- static_cast<RealType>(100), // degrees of freedom
- static_cast<RealType>(16), // non centrality
- static_cast<RealType>(82.35814), // Chi Squared statistic
- static_cast<RealType>(0.1184348822747824e-1), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.1184348822747824e-1), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(300), // degrees of freedom
- static_cast<RealType>(16), // non centrality
- static_cast<RealType>(331.78852), // Chi Squared statistic
- static_cast<RealType>(0.7355956710306709), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.7355956710306709), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(500), // degrees of freedom
- static_cast<RealType>(21), // non centrality
- static_cast<RealType>(459.92612), // Chi Squared statistic
- static_cast<RealType>(0.2797023600800060e-1), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.2797023600800060e-1), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(1), // degrees of freedom
- static_cast<RealType>(1), // non centrality
- static_cast<RealType>(0.00016), // Chi Squared statistic
- static_cast<RealType>(0.6121428929881423e-2), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.6121428929881423e-2), // Q = 1 - P
- tolerance);
- test_spot(
- static_cast<RealType>(1), // degrees of freedom
- static_cast<RealType>(1), // non centrality
- static_cast<RealType>(0.00393), // Chi Squared statistic
- static_cast<RealType>(0.3033814229753780e-1), // Probability of result (CDF), P
- static_cast<RealType>(1 - 0.3033814229753780e-1), // Q = 1 - P
- tolerance);
- RealType tol2 = boost::math::tools::epsilon<RealType>() * 5 * 100; // 5 eps as a percentage
- boost::math::non_central_chi_squared_distribution<RealType> dist(static_cast<RealType>(8), static_cast<RealType>(12));
- RealType x = 7;
- using namespace std; // ADL of std names.
- // mean:
- BOOST_CHECK_CLOSE(
- mean(dist)
- , static_cast<RealType>(8 + 12), tol2);
- // variance:
- BOOST_CHECK_CLOSE(
- variance(dist)
- , static_cast<RealType>(64), tol2);
- // std deviation:
- BOOST_CHECK_CLOSE(
- standard_deviation(dist)
- , static_cast<RealType>(8), tol2);
- // hazard:
- 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);
- // mode:
- BOOST_CHECK_CLOSE(
- mode(dist)
- , static_cast<RealType>(17.184201184730857030170788677340294070728990862663L), sqrt(tolerance * 500));
- BOOST_CHECK_CLOSE(
- median(dist),
- quantile(
- boost::math::non_central_chi_squared_distribution<RealType>(
- static_cast<RealType>(8),
- static_cast<RealType>(12)),
- static_cast<RealType>(0.5)), static_cast<RealType>(tol2));
- // skewness:
- BOOST_CHECK_CLOSE(
- skewness(dist)
- , static_cast<RealType>(0.6875), tol2);
- // kurtosis:
- BOOST_CHECK_CLOSE(
- kurtosis(dist)
- , static_cast<RealType>(3.65625), tol2);
- // kurtosis excess:
- BOOST_CHECK_CLOSE(
- kurtosis_excess(dist)
- , static_cast<RealType>(0.65625), tol2);
- // Error handling checks:
- check_out_of_range<boost::math::non_central_chi_squared_distribution<RealType> >(1, 1);
- BOOST_MATH_CHECK_THROW(pdf(boost::math::non_central_chi_squared_distribution<RealType>(0, 1), 0), std::domain_error);
- BOOST_MATH_CHECK_THROW(pdf(boost::math::non_central_chi_squared_distribution<RealType>(-1, 1), 0), std::domain_error);
- BOOST_MATH_CHECK_THROW(pdf(boost::math::non_central_chi_squared_distribution<RealType>(1, -1), 0), std::domain_error);
- BOOST_MATH_CHECK_THROW(quantile(boost::math::non_central_chi_squared_distribution<RealType>(1, 1), -1), std::domain_error);
- BOOST_MATH_CHECK_THROW(quantile(boost::math::non_central_chi_squared_distribution<RealType>(1, 1), 2), std::domain_error);
- #endif
- } // template <class RealType>void test_spots(RealType)
- template <class T>
- T nccs_cdf(T df, T nc, T x)
- {
- return cdf(boost::math::non_central_chi_squared_distribution<T>(df, nc), x);
- }
- template <class T>
- T nccs_ccdf(T df, T nc, T x)
- {
- return cdf(complement(boost::math::non_central_chi_squared_distribution<T>(df, nc), x));
- }
- template <typename Real, typename T>
- void do_test_nc_chi_squared(T& data, const char* type_name, const char* test)
- {
- typedef Real value_type;
- std::cout << "Testing: " << test << std::endl;
- #ifdef NC_CHI_SQUARED_CDF_FUNCTION_TO_TEST
- value_type(*fp1)(value_type, value_type, value_type) = NC_CHI_SQUARED_CDF_FUNCTION_TO_TEST;
- #else
- value_type(*fp1)(value_type, value_type, value_type) = nccs_cdf;
- #endif
- boost::math::tools::test_result<value_type> result;
- #if !(defined(ERROR_REPORTING_MODE) && !defined(NC_CHI_SQUARED_CDF_FUNCTION_TO_TEST))
- result = boost::math::tools::test_hetero<Real>(
- data,
- bind_func<Real>(fp1, 0, 1, 2),
- extract_result<Real>(3));
- handle_test_result(result, data[result.worst()], result.worst(),
- type_name, "non central chi squared CDF", test);
- #endif
- #if !(defined(ERROR_REPORTING_MODE) && !defined(NC_CHI_SQUARED_CCDF_FUNCTION_TO_TEST))
- #ifdef NC_CHI_SQUARED_CCDF_FUNCTION_TO_TEST
- fp1 = NC_CHI_SQUARED_CCDF_FUNCTION_TO_TEST;
- #else
- fp1 = nccs_ccdf;
- #endif
- result = boost::math::tools::test_hetero<Real>(
- data,
- bind_func<Real>(fp1, 0, 1, 2),
- extract_result<Real>(4));
- handle_test_result(result, data[result.worst()], result.worst(),
- type_name, "non central chi squared CDF complement", test);
- std::cout << std::endl;
- #endif
- }
- template <typename Real, typename T>
- void quantile_sanity_check(T& data, const char* type_name, const char* test)
- {
- #ifndef ERROR_REPORTING_MODE
- typedef Real value_type;
- //
- // Tests with type real_concept take rather too long to run, so
- // for now we'll disable them:
- //
- if(!boost::is_floating_point<value_type>::value)
- return;
- std::cout << "Testing: " << type_name << " quantile sanity check, with tests " << test << std::endl;
- //
- // These sanity checks test for a round trip accuracy of one half
- // of the bits in T, unless T is type float, in which case we check
- // for just one decimal digit. The problem here is the sensitivity
- // of the functions, not their accuracy. This test data was generated
- // for the forward functions, which means that when it is used as
- // the input to the inverses then it is necessarily inexact. This rounding
- // of the input is what makes the data unsuitable for use as an accuracy check,
- // and also demonstrates that you can't in general round-trip these functions.
- // It is however a useful sanity check.
- //
- value_type precision = static_cast<value_type>(ldexp(1.0, 1 - boost::math::policies::digits<value_type, boost::math::policies::policy<> >() / 2)) * 100;
- if(boost::math::policies::digits<value_type, boost::math::policies::policy<> >() < 50)
- precision = 1; // 1% or two decimal digits, all we can hope for when the input is truncated to float
- for(unsigned i = 0; i < data.size(); ++i)
- {
- if(Real(data[i][3]) == 0)
- {
- BOOST_CHECK(0 == quantile(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][3]));
- }
- else if(data[i][3] < 0.9999f)
- {
- value_type p = quantile(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][3]);
- value_type pt = data[i][2];
- BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
- }
- if(data[i][4] == 0)
- {
- BOOST_CHECK(0 == quantile(complement(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][3])));
- }
- else if(data[i][4] < 0.9999f)
- {
- value_type p = quantile(complement(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][4]));
- value_type pt = data[i][2];
- BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
- }
- if(boost::math::tools::digits<value_type>() > 50)
- {
- //
- // Sanity check mode, the accuracy of
- // the mode is at *best* the square root of the accuracy of the PDF:
- //
- #ifndef BOOST_NO_EXCEPTIONS
- try{
- value_type m = mode(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]));
- value_type p = pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m);
- BOOST_CHECK_EX(pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m * (1 + sqrt(precision) * 50)) <= p, i);
- BOOST_CHECK_EX(pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m * (1 - sqrt(precision)) * 50) <= p, i);
- }
- catch(const boost::math::evaluation_error&) {}
- #endif
- //
- // Sanity check degrees-of-freedom finder, don't bother at float
- // precision though as there's not enough data in the probability
- // values to get back to the correct degrees of freedom or
- // non-cenrality parameter:
- //
- #ifndef BOOST_NO_EXCEPTIONS
- try{
- #endif
- if((data[i][3] < 0.99) && (data[i][3] != 0))
- {
- BOOST_CHECK_CLOSE_EX(
- boost::math::non_central_chi_squared_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], data[i][3]),
- data[i][0], precision, i);
- BOOST_CHECK_CLOSE_EX(
- boost::math::non_central_chi_squared_distribution<value_type>::find_non_centrality(data[i][0], data[i][2], data[i][3]),
- data[i][1], precision, i);
- }
- if((data[i][4] < 0.99) && (data[i][4] != 0))
- {
- BOOST_CHECK_CLOSE_EX(
- boost::math::non_central_chi_squared_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], data[i][4])),
- data[i][0], precision, i);
- BOOST_CHECK_CLOSE_EX(
- boost::math::non_central_chi_squared_distribution<value_type>::find_non_centrality(boost::math::complement(data[i][0], data[i][2], data[i][4])),
- data[i][1], precision, i);
- }
- #ifndef BOOST_NO_EXCEPTIONS
- }
- catch(const std::exception& e)
- {
- BOOST_ERROR(e.what());
- }
- #endif
- }
- }
- #endif
- }
- template <typename T>
- void test_accuracy(T, const char* type_name)
- {
- #include "nccs.ipp"
- do_test_nc_chi_squared<T>(nccs, type_name, "Non Central Chi Squared, medium parameters");
- quantile_sanity_check<T>(nccs, type_name, "Non Central Chi Squared, medium parameters");
- #include "nccs_big.ipp"
- do_test_nc_chi_squared<T>(nccs_big, type_name, "Non Central Chi Squared, large parameters");
- quantile_sanity_check<T>(nccs_big, type_name, "Non Central Chi Squared, large parameters");
- }
|