libs/math/test/test_nc_chi_squared.hpp

   1 //  (C) Copyright John Maddock 2007.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
   7 #include <boost/math/concepts/real_concept.hpp>
   8 #define BOOST_TEST_MAIN
   9 #include <boost/test/unit_test.hpp>
  10 #include <boost/test/tools/floating_point_comparison.hpp>
  11 #include <boost/math/distributions/non_central_chi_squared.hpp>
  12 #include <boost/type_traits/is_floating_point.hpp>
  13 #include <boost/array.hpp>
  14 #include "functor.hpp"
  15
  16 #include "handle_test_result.hpp"
  17 #include "table_type.hpp"
  18
  19 #include <iostream>
  20 #include <iomanip>
  21
  22 #define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \
  23       {\
  24       unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
  25       BOOST_CHECK_CLOSE(a, b, prec); \
  26       if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
  27             {\
  28          std::cerr << "Failure was at row " << i << std::endl;\
  29          std::cerr << std::setprecision(35); \
  30          std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
  31          std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
  32             }\
  33       }
  34
  35 #define BOOST_CHECK_EX(a, i) \
  36       {\
  37       unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
  38       BOOST_CHECK(a); \
  39       if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
  40             {\
  41          std::cerr << "Failure was at row " << i << std::endl;\
  42          std::cerr << std::setprecision(35); \
  43          std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
  44          std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
  45             }\
  46       }
  47
  48 template <class RealType>
  49 RealType naive_pdf(RealType v, RealType lam, RealType x)
  50 {
  51    // Formula direct from
  52    // http://mathworld.wolfram.com/NoncentralChi-SquaredDistribution.html
  53    // with no simplification:
  54    RealType sum, term, prefix(1);
  55    RealType eps = boost::math::tools::epsilon<RealType>();
  56    term = sum = pdf(boost::math::chi_squared_distribution<RealType>(v), x);
  57    for(int i = 1;; ++i)
  58    {
  59       prefix *= lam / (2 * i);
  60       term = prefix * pdf(boost::math::chi_squared_distribution<RealType>(v + 2 * i), x);
  61       sum += term;
  62       if(term / sum < eps)
  63          break;
  64    }
  65    return sum * exp(-lam / 2);
  66 }
  67
  68 template <class RealType>
  69 void test_spot(
  70    RealType df,    // Degrees of freedom
  71    RealType ncp,   // non-centrality param
  72    RealType cs,    // Chi Square statistic
  73    RealType P,     // CDF
  74    RealType Q,     // Complement of CDF
  75    RealType tol)   // Test tolerance
  76 {
  77    boost::math::non_central_chi_squared_distribution<RealType> dist(df, ncp);
  78    BOOST_CHECK_CLOSE(
  79       cdf(dist, cs), P, tol);
  80 #ifndef BOOST_NO_EXCEPTIONS
  81    try{
  82       BOOST_CHECK_CLOSE(
  83          pdf(dist, cs), naive_pdf(dist.degrees_of_freedom(), ncp, cs), tol * 150);
  84    }
  85    catch(const std::overflow_error&)
  86    {
  87    }
  88 #endif
  89    if((P < 0.99) && (Q < 0.99))
  90    {
  91       //
  92       // We can only check this if P is not too close to 1,
  93       // so that we can guarantee Q is reasonably free of error:
  94       //
  95       BOOST_CHECK_CLOSE(
  96          cdf(complement(dist, cs)), Q, tol);
  97       BOOST_CHECK_CLOSE(
  98          quantile(dist, P), cs, tol * 10);
  99       BOOST_CHECK_CLOSE(
 100          quantile(complement(dist, Q)), cs, tol * 10);
 101       BOOST_CHECK_CLOSE(
 102          dist.find_degrees_of_freedom(ncp, cs, P), df, tol * 10);
 103       BOOST_CHECK_CLOSE(
 104          dist.find_degrees_of_freedom(boost::math::complement(ncp, cs, Q)), df, tol * 10);
 105       BOOST_CHECK_CLOSE(
 106          dist.find_non_centrality(df, cs, P), ncp, tol * 10);
 107       BOOST_CHECK_CLOSE(
 108          dist.find_non_centrality(boost::math::complement(df, cs, Q)), ncp, tol * 10);
 109    }
 110 }
 111
 112 template <class RealType> // Any floating-point type RealType.
 113 void test_spots(RealType)
 114 {
 115 #ifndef ERROR_REPORTING_MODE
 116    RealType tolerance = (std::max)(
 117       boost::math::tools::epsilon<RealType>(),
 118       (RealType)boost::math::tools::epsilon<double>() * 5) * 150;
 119    //
 120    // At float precision we need to up the tolerance, since
 121    // the input values are rounded off to inexact quantities
 122    // the results get thrown off by a noticeable amount.
 123    //
 124    if(boost::math::tools::digits<RealType>() < 50)
 125       tolerance *= 50;
 126    if(boost::is_floating_point<RealType>::value != 1)
 127       tolerance *= 20; // real_concept special functions are less accurate
 128
 129    std::cout << "Tolerance = " << tolerance << "%." << std::endl;
 130
 131    using boost::math::chi_squared_distribution;
 132    using  ::boost::math::chi_squared;
 133    using  ::boost::math::cdf;
 134    using  ::boost::math::pdf;
 135    //
 136    // Test against the data from Table 6 of:
 137    //
 138    // "Self-Validating Computations of Probabilities for Selected
 139    // Central and Noncentral Univariate Probability Functions."
 140    // Morgan C. Wang; William J. Kennedy
 141    // Journal of the American Statistical Association,
 142    // Vol. 89, No. 427. (Sep., 1994), pp. 878-887.
 143    //
 144    test_spot(
 145       static_cast<RealType>(1),   // degrees of freedom
 146       static_cast<RealType>(6),   // non centrality
 147       static_cast<RealType>(0.00393),   // Chi Squared statistic
 148       static_cast<RealType>(0.2498463724258039e-2),       // Probability of result (CDF), P
 149       static_cast<RealType>(1 - 0.2498463724258039e-2),           // Q = 1 - P
 150       tolerance);
 151    test_spot(
 152       static_cast<RealType>(5),   // degrees of freedom
 153       static_cast<RealType>(1),   // non centrality
 154       static_cast<RealType>(9.23636),   // Chi Squared statistic
 155       static_cast<RealType>(0.8272918751175548),       // Probability of result (CDF), P
 156       static_cast<RealType>(1 - 0.8272918751175548),           // Q = 1 - P
 157       tolerance);
 158    test_spot(
 159       static_cast<RealType>(11),   // degrees of freedom
 160       static_cast<RealType>(21),   // non centrality
 161       static_cast<RealType>(24.72497),   // Chi Squared statistic
 162       static_cast<RealType>(0.2539481822183126),       // Probability of result (CDF), P
 163       static_cast<RealType>(1 - 0.2539481822183126),           // Q = 1 - P
 164       tolerance);
 165    test_spot(
 166       static_cast<RealType>(31),   // degrees of freedom
 167       static_cast<RealType>(6),   // non centrality
 168       static_cast<RealType>(44.98534),   // Chi Squared statistic
 169       static_cast<RealType>(0.8125198785064969),       // Probability of result (CDF), P
 170       static_cast<RealType>(1 - 0.8125198785064969),           // Q = 1 - P
 171       tolerance);
 172    test_spot(
 173       static_cast<RealType>(51),   // degrees of freedom
 174       static_cast<RealType>(1),   // non centrality
 175       static_cast<RealType>(38.56038),   // Chi Squared statistic
 176       static_cast<RealType>(0.8519497361859118e-1),       // Probability of result (CDF), P
 177       static_cast<RealType>(1 - 0.8519497361859118e-1),           // Q = 1 - P
 178       tolerance * 2);
 179    test_spot(
 180       static_cast<RealType>(100),   // degrees of freedom
 181       static_cast<RealType>(16),   // non centrality
 182       static_cast<RealType>(82.35814),   // Chi Squared statistic
 183       static_cast<RealType>(0.1184348822747824e-1),       // Probability of result (CDF), P
 184       static_cast<RealType>(1 - 0.1184348822747824e-1),           // Q = 1 - P
 185       tolerance);
 186    test_spot(
 187       static_cast<RealType>(300),   // degrees of freedom
 188       static_cast<RealType>(16),   // non centrality
 189       static_cast<RealType>(331.78852),   // Chi Squared statistic
 190       static_cast<RealType>(0.7355956710306709),       // Probability of result (CDF), P
 191       static_cast<RealType>(1 - 0.7355956710306709),           // Q = 1 - P
 192       tolerance);
 193    test_spot(
 194       static_cast<RealType>(500),   // degrees of freedom
 195       static_cast<RealType>(21),   // non centrality
 196       static_cast<RealType>(459.92612),   // Chi Squared statistic
 197       static_cast<RealType>(0.2797023600800060e-1),       // Probability of result (CDF), P
 198       static_cast<RealType>(1 - 0.2797023600800060e-1),           // Q = 1 - P
 199       tolerance);
 200    test_spot(
 201       static_cast<RealType>(1),   // degrees of freedom
 202       static_cast<RealType>(1),   // non centrality
 203       static_cast<RealType>(0.00016),   // Chi Squared statistic
 204       static_cast<RealType>(0.6121428929881423e-2),       // Probability of result (CDF), P
 205       static_cast<RealType>(1 - 0.6121428929881423e-2),           // Q = 1 - P
 206       tolerance);
 207    test_spot(
 208       static_cast<RealType>(1),   // degrees of freedom
 209       static_cast<RealType>(1),   // non centrality
 210       static_cast<RealType>(0.00393),   // Chi Squared statistic
 211       static_cast<RealType>(0.3033814229753780e-1),       // Probability of result (CDF), P
 212       static_cast<RealType>(1 - 0.3033814229753780e-1),           // Q = 1 - P
 213       tolerance);
 214
 215    RealType tol2 = boost::math::tools::epsilon<RealType>() * 5 * 100; // 5 eps as a percentage
 216    boost::math::non_central_chi_squared_distribution<RealType> dist(static_cast<RealType>(8), static_cast<RealType>(12));
 217    RealType x = 7;
 218    using namespace std; // ADL of std names.
 219    // mean:
 220    BOOST_CHECK_CLOSE(
 221       mean(dist)
 222       , static_cast<RealType>(8 + 12), tol2);
 223    // variance:
 224    BOOST_CHECK_CLOSE(
 225       variance(dist)
 226       , static_cast<RealType>(64), tol2);
 227    // std deviation:
 228    BOOST_CHECK_CLOSE(
 229       standard_deviation(dist)
 230       , static_cast<RealType>(8), tol2);
 231    // hazard:
 232    BOOST_CHECK_CLOSE(
 233       hazard(dist, x)
 234       , pdf(dist, x) / cdf(complement(dist, x)), tol2);
 235    // cumulative hazard:
 236    BOOST_CHECK_CLOSE(
 237       chf(dist, x)
 238       , -log(cdf(complement(dist, x))), tol2);
 239    // coefficient_of_variation:
 240    BOOST_CHECK_CLOSE(
 241       coefficient_of_variation(dist)
 242       , standard_deviation(dist) / mean(dist), tol2);
 243    // mode:
 244    BOOST_CHECK_CLOSE(
 245       mode(dist)
 246       , static_cast<RealType>(17.184201184730857030170788677340294070728990862663L), sqrt(tolerance * 500));
 247    BOOST_CHECK_CLOSE(
 248       median(dist),
 249       quantile(
 250       boost::math::non_central_chi_squared_distribution<RealType>(
 251       static_cast<RealType>(8),
 252       static_cast<RealType>(12)),
 253       static_cast<RealType>(0.5)), static_cast<RealType>(tol2));
 254    // skewness:
 255    BOOST_CHECK_CLOSE(
 256       skewness(dist)
 257       , static_cast<RealType>(0.6875), tol2);
 258    // kurtosis:
 259    BOOST_CHECK_CLOSE(
 260       kurtosis(dist)
 261       , static_cast<RealType>(3.65625), tol2);
 262    // kurtosis excess:
 263    BOOST_CHECK_CLOSE(
 264       kurtosis_excess(dist)
 265       , static_cast<RealType>(0.65625), tol2);
 266
 267    // Error handling checks:
 268    check_out_of_range<boost::math::non_central_chi_squared_distribution<RealType> >(1, 1);
 269    BOOST_MATH_CHECK_THROW(pdf(boost::math::non_central_chi_squared_distribution<RealType>(0, 1), 0), std::domain_error);
 270    BOOST_MATH_CHECK_THROW(pdf(boost::math::non_central_chi_squared_distribution<RealType>(-1, 1), 0), std::domain_error);
 271    BOOST_MATH_CHECK_THROW(pdf(boost::math::non_central_chi_squared_distribution<RealType>(1, -1), 0), std::domain_error);
 272    BOOST_MATH_CHECK_THROW(quantile(boost::math::non_central_chi_squared_distribution<RealType>(1, 1), -1), std::domain_error);
 273    BOOST_MATH_CHECK_THROW(quantile(boost::math::non_central_chi_squared_distribution<RealType>(1, 1), 2), std::domain_error);
 274 #endif
 275 } // template <class RealType>void test_spots(RealType)
 276
 277 template <class T>
 278 T nccs_cdf(T df, T nc, T x)
 279 {
 280    return cdf(boost::math::non_central_chi_squared_distribution<T>(df, nc), x);
 281 }
 282
 283 template <class T>
 284 T nccs_ccdf(T df, T nc, T x)
 285 {
 286    return cdf(complement(boost::math::non_central_chi_squared_distribution<T>(df, nc), x));
 287 }
 288
 289 template <typename Real, typename T>
 290 void do_test_nc_chi_squared(T& data, const char* type_name, const char* test)
 291 {
 292    typedef Real                   value_type;
 293
 294    std::cout << "Testing: " << test << std::endl;
 295
 296 #ifdef NC_CHI_SQUARED_CDF_FUNCTION_TO_TEST
 297    value_type(*fp1)(value_type, value_type, value_type) = NC_CHI_SQUARED_CDF_FUNCTION_TO_TEST;
 298 #else
 299    value_type(*fp1)(value_type, value_type, value_type) = nccs_cdf;
 300 #endif
 301    boost::math::tools::test_result<value_type> result;
 302
 303 #if !(defined(ERROR_REPORTING_MODE) && !defined(NC_CHI_SQUARED_CDF_FUNCTION_TO_TEST))
 304    result = boost::math::tools::test_hetero<Real>(
 305       data,
 306       bind_func<Real>(fp1, 0, 1, 2),
 307       extract_result<Real>(3));
 308    handle_test_result(result, data[result.worst()], result.worst(),
 309       type_name, "non central chi squared CDF", test);
 310 #endif
 311 #if !(defined(ERROR_REPORTING_MODE) && !defined(NC_CHI_SQUARED_CCDF_FUNCTION_TO_TEST))
 312 #ifdef NC_CHI_SQUARED_CCDF_FUNCTION_TO_TEST
 313    fp1 = NC_CHI_SQUARED_CCDF_FUNCTION_TO_TEST;
 314 #else
 315    fp1 = nccs_ccdf;
 316 #endif
 317    result = boost::math::tools::test_hetero<Real>(
 318       data,
 319       bind_func<Real>(fp1, 0, 1, 2),
 320       extract_result<Real>(4));
 321    handle_test_result(result, data[result.worst()], result.worst(),
 322       type_name, "non central chi squared CDF complement", test);
 323
 324    std::cout << std::endl;
 325 #endif
 326 }
 327
 328 template <typename Real, typename T>
 329 void quantile_sanity_check(T& data, const char* type_name, const char* test)
 330 {
 331 #ifndef ERROR_REPORTING_MODE
 332    typedef Real                   value_type;
 333
 334    //
 335    // Tests with type real_concept take rather too long to run, so
 336    // for now we'll disable them:
 337    //
 338    if(!boost::is_floating_point<value_type>::value)
 339       return;
 340
 341    std::cout << "Testing: " << type_name << " quantile sanity check, with tests " << test << std::endl;
 342
 343    //
 344    // These sanity checks test for a round trip accuracy of one half
 345    // of the bits in T, unless T is type float, in which case we check
 346    // for just one decimal digit.  The problem here is the sensitivity
 347    // of the functions, not their accuracy.  This test data was generated
 348    // for the forward functions, which means that when it is used as
 349    // the input to the inverses then it is necessarily inexact.  This rounding
 350    // of the input is what makes the data unsuitable for use as an accuracy check,
 351    // and also demonstrates that you can't in general round-trip these functions.
 352    // It is however a useful sanity check.
 353    //
 354    value_type precision = static_cast<value_type>(ldexp(1.0, 1 - boost::math::policies::digits<value_type, boost::math::policies::policy<> >() / 2)) * 100;
 355    if(boost::math::policies::digits<value_type, boost::math::policies::policy<> >() < 50)
 356       precision = 1;   // 1% or two decimal digits, all we can hope for when the input is truncated to float
 357
 358    for(unsigned i = 0; i < data.size(); ++i)
 359    {
 360       if(Real(data[i][3]) == 0)
 361       {
 362          BOOST_CHECK(0 == quantile(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][3]));
 363       }
 364       else if(data[i][3] < 0.9999f)
 365       {
 366          value_type p = quantile(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][3]);
 367          value_type pt = data[i][2];
 368          BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
 369       }
 370       if(data[i][4] == 0)
 371       {
 372          BOOST_CHECK(0 == quantile(complement(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][3])));
 373       }
 374       else if(data[i][4] < 0.9999f)
 375       {
 376          value_type p = quantile(complement(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][4]));
 377          value_type pt = data[i][2];
 378          BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
 379       }
 380       if(boost::math::tools::digits<value_type>() > 50)
 381       {
 382          //
 383          // Sanity check mode, the accuracy of
 384          // the mode is at *best* the square root of the accuracy of the PDF:
 385          //
 386 #ifndef BOOST_NO_EXCEPTIONS
 387          try{
 388             value_type m = mode(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]));
 389             value_type p = pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m);
 390             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);
 391             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);
 392          }
 393          catch(const boost::math::evaluation_error&) {}
 394 #endif
 395          //
 396          // Sanity check degrees-of-freedom finder, don't bother at float
 397          // precision though as there's not enough data in the probability
 398          // values to get back to the correct degrees of freedom or
 399          // non-cenrality parameter:
 400          //
 401 #ifndef BOOST_NO_EXCEPTIONS
 402          try{
 403 #endif
 404             if((data[i][3] < 0.99) && (data[i][3] != 0))
 405             {
 406                BOOST_CHECK_CLOSE_EX(
 407                   boost::math::non_central_chi_squared_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], data[i][3]),
 408                   data[i][0], precision, i);
 409                BOOST_CHECK_CLOSE_EX(
 410                   boost::math::non_central_chi_squared_distribution<value_type>::find_non_centrality(data[i][0], data[i][2], data[i][3]),
 411                   data[i][1], precision, i);
 412             }
 413             if((data[i][4] < 0.99) && (data[i][4] != 0))
 414             {
 415                BOOST_CHECK_CLOSE_EX(
 416                   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])),
 417                   data[i][0], precision, i);
 418                BOOST_CHECK_CLOSE_EX(
 419                   boost::math::non_central_chi_squared_distribution<value_type>::find_non_centrality(boost::math::complement(data[i][0], data[i][2], data[i][4])),
 420                   data[i][1], precision, i);
 421             }
 422 #ifndef BOOST_NO_EXCEPTIONS
 423          }
 424          catch(const std::exception& e)
 425          {
 426             BOOST_ERROR(e.what());
 427          }
 428 #endif
 429       }
 430    }
 431 #endif
 432 }
 433
 434 template <typename T>
 435 void test_accuracy(T, const char* type_name)
 436 {
 437 #include "nccs.ipp"
 438    do_test_nc_chi_squared<T>(nccs, type_name, "Non Central Chi Squared, medium parameters");
 439    quantile_sanity_check<T>(nccs, type_name, "Non Central Chi Squared, medium parameters");
 440
 441 #include "nccs_big.ipp"
 442    do_test_nc_chi_squared<T>(nccs_big, type_name, "Non Central Chi Squared, large parameters");
 443    quantile_sanity_check<T>(nccs_big, type_name, "Non Central Chi Squared, large parameters");
 444 }