| // test_nc_t.cpp |
| |
| // Copyright John Maddock 2008. |
| |
| // 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) |
| |
| #include <pch.hpp> |
| |
| #ifdef _MSC_VER |
| #pragma warning (disable:4127 4512) |
| #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/concepts/real_concept.hpp> // for real_concept |
| #include <boost/math/distributions/non_central_t.hpp> // for chi_squared_distribution |
| #include <boost/test/test_exec_monitor.hpp> // for test_main |
| #include <boost/test/results_collector.hpp> |
| #include <boost/test/unit_test.hpp> |
| #include <boost/test/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE |
| |
| #include "functor.hpp" |
| #include "handle_test_result.hpp" |
| |
| #include <iostream> |
| using std::cout; |
| using std::endl; |
| #include <limits> |
| using std::numeric_limits; |
| |
| #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;\ |
| }\ |
| } |
| |
| void expected_results() |
| { |
| // |
| // Define the max and mean errors expected for |
| // various compilers and platforms. |
| // |
| const char* largest_type; |
| #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
| if(boost::math::policies::digits<double, boost::math::policies::policy<> >() == boost::math::policies::digits<long double, boost::math::policies::policy<> >()) |
| { |
| largest_type = "(long\\s+)?double|real_concept"; |
| } |
| else |
| { |
| largest_type = "long double|real_concept"; |
| } |
| #else |
| largest_type = "(long\\s+)?double|real_concept"; |
| #endif |
| |
| // |
| // Catch all cases come last: |
| // |
| add_expected_result( |
| "[^|]*", // compiler |
| "[^|]*", // stdlib |
| "[^|]*", // platform |
| "real_concept", // test type(s) |
| "[^|]*", // test data group |
| "[^|]*", 300000, 100000); // test function |
| add_expected_result( |
| "[^|]*", // compiler |
| "[^|]*", // stdlib |
| "[^|]*", // platform |
| largest_type, // test type(s) |
| "[^|]*", // test data group |
| "[^|]*", 250, 50); // test function |
| |
| // |
| // Finish off by printing out the compiler/stdlib/platform names, |
| // we do this to make it easier to mark up expected error rates. |
| // |
| std::cout << "Tests run with " << BOOST_COMPILER << ", " |
| << BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl; |
| } |
| |
| template <class RealType> |
| RealType naive_pdf(RealType v, RealType delta, RealType x) |
| { |
| } |
| |
| template <class RealType> |
| RealType naive_mean(RealType v, RealType delta) |
| { |
| using boost::math::tgamma; |
| return delta * sqrt(v / 2) * tgamma((v-1)/2) / tgamma(v/2); |
| } |
| |
| float naive_mean(float v, float delta) |
| { |
| return (float)naive_mean((double)v, (double)delta); |
| } |
| |
| template <class RealType> |
| RealType naive_variance(RealType v, RealType delta) |
| { |
| using boost::math::tgamma; |
| RealType r = tgamma((v-1)/2) / tgamma(v/2); |
| r *= r; |
| r *= -delta * delta * v / 2; |
| r += (1 + delta * delta) * v / (v - 2); |
| return r; |
| } |
| |
| float naive_variance(float v, float delta) |
| { |
| return (float)naive_variance((double)v, (double)delta); |
| } |
| |
| template <class RealType> |
| RealType naive_skewness(RealType v, RealType delta) |
| { |
| using boost::math::tgamma; |
| RealType tgr = tgamma((v-1)/2) / tgamma(v / 2); |
| RealType r = delta * sqrt(v) * tgamma((v-1)/2) |
| * (v * (-3 + delta * delta + 2 * v) / ((-3 + v) * (-2 + v)) |
| - 2 * ((1 + delta * delta) * v / (-2 + v) - delta * delta * v * tgr * tgr / 2)); |
| r /= boost::math::constants::root_two<RealType>() |
| * pow(((1+delta*delta) * v / (-2+v) - delta*delta*v*tgr*tgr/2), RealType(1.5f)) |
| * tgamma(v/2); |
| return r; |
| } |
| |
| float naive_skewness(float v, float delta) |
| { |
| return (float)naive_skewness((double)v, (double)delta); |
| } |
| |
| template <class RealType> |
| RealType naive_kurtosis_excess(RealType v, RealType delta) |
| { |
| using boost::math::tgamma; |
| RealType tgr = tgamma((v-1)/2) / tgamma(v / 2); |
| RealType r = -delta * delta * v * tgr * tgr / 2; |
| r *= v * (delta * delta * (1 + v) + 3 * (-5 + 3 * v)) / ((-3 + v)*(-2+v)) |
| - 3 * ((1 + delta * delta) * v / (-2 + v) - delta * delta * v * tgr * tgr / 2); |
| r += (3 + 6 * delta * delta + delta * delta * delta * delta)* v * v |
| / ((-4+v) * (-2+v)); |
| r /= (1+delta*delta)*v / (-2+v) - delta*delta*v *tgr*tgr/2; |
| r /= (1+delta*delta)*v / (-2+v) - delta*delta*v *tgr*tgr/2; |
| return r; |
| } |
| |
| float naive_kurtosis_excess(float v, float delta) |
| { |
| return (float)naive_kurtosis_excess((double)v, (double)delta); |
| } |
| |
| template <class RealType> |
| void test_spot( |
| RealType df, // Degrees of freedom |
| RealType ncp, // non-centrality param |
| RealType t, // T statistic |
| RealType P, // CDF |
| RealType Q, // Complement of CDF |
| RealType tol) // Test tolerance |
| { |
| boost::math::non_central_t_distribution<RealType> dist(df, ncp); |
| BOOST_CHECK_CLOSE( |
| cdf(dist, t), P, tol); |
| try{ |
| BOOST_CHECK_CLOSE( |
| mean(dist), naive_mean(df, ncp), tol); |
| BOOST_CHECK_CLOSE( |
| variance(dist), naive_variance(df, ncp), tol); |
| BOOST_CHECK_CLOSE( |
| skewness(dist), naive_skewness(df, ncp), tol * 10); |
| BOOST_CHECK_CLOSE( |
| kurtosis_excess(dist), naive_kurtosis_excess(df, ncp), tol * 50); |
| BOOST_CHECK_CLOSE( |
| kurtosis(dist), 3 + naive_kurtosis_excess(df, ncp), tol * 50); |
| } |
| catch(const std::domain_error&) |
| { |
| } |
| /* |
| BOOST_CHECK_CLOSE( |
| pdf(dist, t), naive_pdf(dist.degrees_of_freedom(), ncp, t), tol * 50); |
| */ |
| if((P < 0.99) && (Q < 0.99)) |
| { |
| // |
| // We can only check this if P is not too close to 1, |
| // so that we can guarentee Q is reasonably free of error: |
| // |
| BOOST_CHECK_CLOSE( |
| cdf(complement(dist, t)), Q, tol); |
| BOOST_CHECK_CLOSE( |
| quantile(dist, P), t, tol * 10); |
| BOOST_CHECK_CLOSE( |
| quantile(complement(dist, Q)), t, tol * 10); |
| /* |
| BOOST_CHECK_CLOSE( |
| dist.find_degrees_of_freedom(ncp, t, P), df, tol * 10); |
| BOOST_CHECK_CLOSE( |
| dist.find_degrees_of_freedom(boost::math::complement(ncp, t, Q)), df, tol * 10); |
| BOOST_CHECK_CLOSE( |
| dist.find_non_centrality(df, t, P), ncp, tol * 10); |
| BOOST_CHECK_CLOSE( |
| dist.find_non_centrality(boost::math::complement(df, t, Q)), ncp, tol * 10); |
| */ |
| } |
| } |
| |
| template <class RealType> // Any floating-point type RealType. |
| void test_spots(RealType) |
| { |
| // |
| // Approx limit of test data is 12 digits expressed here as a persentage: |
| // |
| RealType tolerance = (std::max)( |
| boost::math::tools::epsilon<RealType>(), |
| (RealType)5e-12f) * 100; |
| // |
| // 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 |
| |
| cout << "Tolerance = " << tolerance << "%." << endl; |
| |
| // |
| // Test data is taken from: |
| // |
| // Computing discrete mixtures of continuous |
| // distributions: noncentral chisquare, noncentral t |
| // and the distribution of the square of the sample |
| // multiple correlation coeficient. |
| // Denise Benton, K. Krishnamoorthy. |
| // Computational Statistics & Data Analysis 43 (2003) 249 - 267 |
| // |
| test_spot( |
| static_cast<RealType>(3), // degrees of freedom |
| static_cast<RealType>(1), // non centrality |
| static_cast<RealType>(2.34), // T |
| static_cast<RealType>(0.801888999613917), // Probability of result (CDF), P |
| static_cast<RealType>(1-0.801888999613917), // Q = 1 - P |
| tolerance); |
| test_spot( |
| static_cast<RealType>(126), // degrees of freedom |
| static_cast<RealType>(-2), // non centrality |
| static_cast<RealType>(-4.33), // T |
| static_cast<RealType>(1.252846196792878e-2), // Probability of result (CDF), P |
| static_cast<RealType>(1-1.252846196792878e-2), // Q = 1 - P |
| tolerance); |
| test_spot( |
| static_cast<RealType>(20), // degrees of freedom |
| static_cast<RealType>(23), // non centrality |
| static_cast<RealType>(23), // T |
| static_cast<RealType>(0.460134400391924), // Probability of result (CDF), P |
| static_cast<RealType>(1-0.460134400391924), // Q = 1 - P |
| tolerance); |
| test_spot( |
| static_cast<RealType>(20), // degrees of freedom |
| static_cast<RealType>(33), // non centrality |
| static_cast<RealType>(34), // T |
| static_cast<RealType>(0.532008386378725), // Probability of result (CDF), P |
| static_cast<RealType>(1-0.532008386378725), // Q = 1 - P |
| tolerance); |
| test_spot( |
| static_cast<RealType>(12), // degrees of freedom |
| static_cast<RealType>(38), // non centrality |
| static_cast<RealType>(39), // T |
| static_cast<RealType>(0.495868184917805), // Probability of result (CDF), P |
| static_cast<RealType>(1-0.495868184917805), // Q = 1 - P |
| tolerance); |
| test_spot( |
| static_cast<RealType>(12), // degrees of freedom |
| static_cast<RealType>(39), // non centrality |
| static_cast<RealType>(39), // T |
| static_cast<RealType>(0.446304024668836), // Probability of result (CDF), P |
| static_cast<RealType>(1-0.446304024668836), // Q = 1 - P |
| tolerance); |
| test_spot( |
| static_cast<RealType>(200), // degrees of freedom |
| static_cast<RealType>(38), // non centrality |
| static_cast<RealType>(39), // T |
| static_cast<RealType>(0.666194209961795), // Probability of result (CDF), P |
| static_cast<RealType>(1-0.666194209961795), // Q = 1 - P |
| tolerance); |
| test_spot( |
| static_cast<RealType>(200), // degrees of freedom |
| static_cast<RealType>(42), // non centrality |
| static_cast<RealType>(40), // T |
| static_cast<RealType>(0.179292265426085), // Probability of result (CDF), P |
| static_cast<RealType>(1-0.179292265426085), // Q = 1 - P |
| tolerance); |
| |
| boost::math::non_central_t_distribution<RealType> dist(static_cast<RealType>(8), static_cast<RealType>(12)); |
| BOOST_CHECK_CLOSE(pdf(dist, 12), static_cast<RealType>(1.235329715425894935157684607751972713457e-1L), tolerance); |
| BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution<RealType>(126, -2), -4), static_cast<RealType>(5.797932289365814702402873546466798025787e-2L), tolerance); |
| BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution<RealType>(126, 2), 4), static_cast<RealType>(5.797932289365814702402873546466798025787e-2L), tolerance); |
| BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution<RealType>(126, 2), 0), static_cast<RealType>(5.388394890639957139696546086044839573749e-2L), tolerance); |
| } // template <class RealType>void test_spots(RealType) |
| |
| template <class T> |
| T nct_cdf(T df, T nc, T x) |
| { |
| return cdf(boost::math::non_central_t_distribution<T>(df, nc), x); |
| } |
| |
| template <class T> |
| T nct_ccdf(T df, T nc, T x) |
| { |
| return cdf(complement(boost::math::non_central_t_distribution<T>(df, nc), x)); |
| } |
| |
| template <typename T> |
| void do_test_nc_t(T& data, const char* type_name, const char* test) |
| { |
| typedef typename T::value_type row_type; |
| typedef typename row_type::value_type value_type; |
| |
| std::cout << "Testing: " << test << std::endl; |
| |
| value_type (*fp1)(value_type, value_type, value_type) = nct_cdf; |
| boost::math::tools::test_result<value_type> result; |
| |
| result = boost::math::tools::test( |
| data, |
| bind_func(fp1, 0, 1, 2), |
| extract_result(3)); |
| handle_test_result(result, data[result.worst()], result.worst(), |
| type_name, "CDF", test); |
| |
| fp1 = nct_ccdf; |
| result = boost::math::tools::test( |
| data, |
| bind_func(fp1, 0, 1, 2), |
| extract_result(4)); |
| handle_test_result(result, data[result.worst()], result.worst(), |
| type_name, "CCDF", test); |
| |
| std::cout << std::endl; |
| |
| } |
| |
| template <typename T> |
| void quantile_sanity_check(T& data, const char* type_name, const char* test) |
| { |
| typedef typename T::value_type row_type; |
| typedef typename row_type::value_type 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(data[i][3] == 0) |
| { |
| BOOST_CHECK(0 == quantile(boost::math::non_central_t_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_t_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_t_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_t_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: |
| // |
| try{ |
| value_type m = mode(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1])); |
| value_type p = pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m); |
| BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m * (1 + sqrt(precision) * 100)) <= p, i); |
| BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m * (1 - sqrt(precision)) * 100) <= p, i); |
| } |
| catch(const boost::math::evaluation_error& ) {} |
| #if 0 |
| // |
| // 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: |
| // |
| try{ |
| if((data[i][3] < 0.99) && (data[i][3] != 0)) |
| { |
| BOOST_CHECK_CLOSE_EX( |
| boost::math::non_central_t_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_t_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_t_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_t_distribution<value_type>::find_non_centrality(boost::math::complement(data[i][0], data[i][2], data[i][4])), |
| data[i][1], precision, i); |
| } |
| } |
| catch(const std::exception& e) |
| { |
| BOOST_ERROR(e.what()); |
| } |
| #endif |
| } |
| } |
| } |
| |
| template <typename T> |
| void test_accuracy(T, const char* type_name) |
| { |
| #include "nct.ipp" |
| do_test_nc_t(nct, type_name, "Non Central T"); |
| quantile_sanity_check(nct, type_name, "Non Central T"); |
| } |
| |
| int test_main(int, char* []) |
| { |
| BOOST_MATH_CONTROL_FP; |
| // Basic sanity-check spot values. |
| expected_results(); |
| |
| // (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 |
| #endif |
| |
| #ifdef TEST_FLOAT |
| test_accuracy(0.0F, "float"); // Test float. |
| #endif |
| #ifdef TEST_DOUBLE |
| test_accuracy(0.0, "double"); // Test double. |
| #endif |
| #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
| #ifdef TEST_LDOUBLE |
| test_accuracy(0.0L, "long double"); // Test long double. |
| #endif |
| #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) |
| #ifdef TEST_REAL_CONCEPT |
| test_accuracy(boost::math::concepts::real_concept(0.), "real_concept"); // Test real concept. |
| #endif |
| #endif |
| #endif |
| return 0; |
| } // int test_main(int, char* []) |
| |
