libs/math/test/test_nc_t.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_t.hpp>
  12 #include <boost/type_traits/is_floating_point.hpp>
  13 #include <boost/array.hpp>
  14 #include "functor.hpp"
  15 #include "test_out_of_range.hpp"
  16
  17 #include "handle_test_result.hpp"
  18 #include "table_type.hpp"
  19
  20 #define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \
  21       {\
  22       unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
  23       BOOST_CHECK_CLOSE(a, b, prec); \
  24       if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
  25             {\
  26          std::cerr << "Failure was at row " << i << std::endl;\
  27          std::cerr << std::setprecision(35); \
  28          std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
  29          std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
  30             }\
  31       }
  32
  33 #define BOOST_CHECK_EX(a, i) \
  34       {\
  35       unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
  36       BOOST_CHECK(a); \
  37       if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
  38             {\
  39          std::cerr << "Failure was at row " << i << std::endl;\
  40          std::cerr << std::setprecision(35); \
  41          std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
  42          std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
  43             }\
  44       }
  45
  46 template <class RealType>
  47 RealType naive_pdf(RealType v, RealType delta, RealType x)
  48 {
  49 }
  50
  51 template <class RealType>
  52 RealType naive_mean(RealType v, RealType delta)
  53 {
  54    using boost::math::tgamma;
  55    return delta * sqrt(v / 2) * tgamma((v - 1) / 2) / tgamma(v / 2);
  56 }
  57
  58 float naive_mean(float v, float delta)
  59 {
  60    return (float)naive_mean((double)v, (double)delta);
  61 }
  62
  63 template <class RealType>
  64 RealType naive_variance(RealType v, RealType delta)
  65 {
  66    using boost::math::tgamma;
  67    RealType r = tgamma((v - 1) / 2) / tgamma(v / 2);
  68    r *= r;
  69    r *= -delta * delta * v / 2;
  70    r += (1 + delta * delta) * v / (v - 2);
  71    return r;
  72 }
  73
  74 float naive_variance(float v, float delta)
  75 {
  76    return (float)naive_variance((double)v, (double)delta);
  77 }
  78
  79 template <class RealType>
  80 RealType naive_skewness(RealType v, RealType delta)
  81 {
  82    using boost::math::tgamma;
  83    RealType tgr = tgamma((v - 1) / 2) / tgamma(v / 2);
  84    RealType r = delta * sqrt(v) * tgamma((v - 1) / 2)
  85       * (v * (-3 + delta * delta + 2 * v) / ((-3 + v) * (-2 + v))
  86       - 2 * ((1 + delta * delta) * v / (-2 + v) - delta * delta * v * tgr * tgr / 2));
  87    r /= boost::math::constants::root_two<RealType>()
  88       * pow(((1 + delta*delta) * v / (-2 + v) - delta*delta*v*tgr*tgr / 2), RealType(1.5f))
  89       * tgamma(v / 2);
  90    return r;
  91 }
  92
  93 float naive_skewness(float v, float delta)
  94 {
  95    return (float)naive_skewness((double)v, (double)delta);
  96 }
  97
  98 template <class RealType>
  99 RealType naive_kurtosis_excess(RealType v, RealType delta)
 100 {
 101    using boost::math::tgamma;
 102    RealType tgr = tgamma((v - 1) / 2) / tgamma(v / 2);
 103    RealType r = -delta * delta * v * tgr * tgr / 2;
 104    r *= v * (delta * delta * (1 + v) + 3 * (-5 + 3 * v)) / ((-3 + v)*(-2 + v))
 105       - 3 * ((1 + delta * delta) * v / (-2 + v) - delta * delta * v * tgr * tgr / 2);
 106    r += (3 + 6 * delta * delta + delta * delta * delta * delta)* v * v
 107       / ((-4 + v) * (-2 + v));
 108    r /= (1 + delta*delta)*v / (-2 + v) - delta*delta*v *tgr*tgr / 2;
 109    r /= (1 + delta*delta)*v / (-2 + v) - delta*delta*v *tgr*tgr / 2;
 110    return r;
 111 }
 112
 113 float naive_kurtosis_excess(float v, float delta)
 114 {
 115    return (float)naive_kurtosis_excess((double)v, (double)delta);
 116 }
 117
 118 template <class RealType>
 119 void test_spot(
 120    RealType df,    // Degrees of freedom
 121    RealType ncp,   // non-centrality param
 122    RealType t,     // T statistic
 123    RealType P,     // CDF
 124    RealType Q,     // Complement of CDF
 125    RealType tol)   // Test tolerance
 126 {
 127    // An extra fudge factor for real_concept which has a less accurate tgamma:
 128    RealType tolerance_tgamma_extra = std::numeric_limits<RealType>::is_specialized ? 1 : 5;
 129
 130    boost::math::non_central_t_distribution<RealType> dist(df, ncp);
 131    BOOST_CHECK_CLOSE(
 132       cdf(dist, t), P, tol);
 133 #ifndef BOOST_NO_EXCEPTIONS
 134    try{
 135       BOOST_CHECK_CLOSE(
 136          mean(dist), naive_mean(df, ncp), tol);
 137       BOOST_CHECK_CLOSE(
 138          variance(dist), naive_variance(df, ncp), tol);
 139       BOOST_CHECK_CLOSE(
 140          skewness(dist), naive_skewness(df, ncp), tol * 10 * tolerance_tgamma_extra);
 141       BOOST_CHECK_CLOSE(
 142          kurtosis_excess(dist), naive_kurtosis_excess(df, ncp), tol * 50 * tolerance_tgamma_extra);
 143       BOOST_CHECK_CLOSE(
 144          kurtosis(dist), 3 + naive_kurtosis_excess(df, ncp), tol * 50 * tolerance_tgamma_extra);
 145    }
 146    catch(const std::domain_error&)
 147    {
 148    }
 149 #endif
 150    /*
 151    BOOST_CHECK_CLOSE(
 152    pdf(dist, t), naive_pdf(dist.degrees_of_freedom(), ncp, t), tol * 50);
 153    */
 154    if((P < 0.99) && (Q < 0.99))
 155    {
 156       //
 157       // We can only check this if P is not too close to 1,
 158       // so that we can guarantee Q is reasonably free of error:
 159       //
 160       BOOST_CHECK_CLOSE(
 161          cdf(complement(dist, t)), Q, tol);
 162       BOOST_CHECK_CLOSE(
 163          quantile(dist, P), t, tol * 10);
 164       BOOST_CHECK_CLOSE(
 165          quantile(complement(dist, Q)), t, tol * 10);
 166       /*  Removed because can give more than one solution.
 167       BOOST_CHECK_CLOSE(
 168       dist.find_degrees_of_freedom(ncp, t, P), df, tol * 10);
 169       BOOST_CHECK_CLOSE(
 170       dist.find_degrees_of_freedom(boost::math::complement(ncp, t, Q)), df, tol * 10);
 171       BOOST_CHECK_CLOSE(
 172       dist.find_non_centrality(df, t, P), ncp, tol * 10);
 173       BOOST_CHECK_CLOSE(
 174       dist.find_non_centrality(boost::math::complement(df, t, Q)), ncp, tol * 10);
 175       */
 176    }
 177 }
 178
 179 template <class RealType> // Any floating-point type RealType.
 180 void test_spots(RealType)
 181 {
 182    using namespace std;
 183    //
 184    // Approx limit of test data is 12 digits expressed here as a percentage:
 185    //
 186    RealType tolerance = (std::max)(
 187       boost::math::tools::epsilon<RealType>(),
 188       (RealType)5e-12f) * 100;
 189    //
 190    // At float precision we need to up the tolerance, since
 191    // the input values are rounded off to inexact quantities
 192    // the results get thrown off by a noticeable amount.
 193    //
 194    if(boost::math::tools::digits<RealType>() < 50)
 195       tolerance *= 50;
 196    if(boost::is_floating_point<RealType>::value != 1)
 197       tolerance *= 20; // real_concept special functions are less accurate
 198
 199    cout << "Tolerance = " << tolerance << "%." << endl;
 200
 201    //
 202    // Test data is taken from:
 203    //
 204    // Computing discrete mixtures of continuous
 205    // distributions: noncentral chisquare, noncentral t
 206    // and the distribution of the square of the sample
 207    // multiple correlation coeficient.
 208    // Denise Benton, K. Krishnamoorthy.
 209    // Computational Statistics & Data Analysis 43 (2003) 249 - 267
 210    //
 211    test_spot(
 212       static_cast<RealType>(3),   // degrees of freedom
 213       static_cast<RealType>(1),   // non centrality
 214       static_cast<RealType>(2.34),   // T
 215       static_cast<RealType>(0.801888999613917),       // Probability of result (CDF), P
 216       static_cast<RealType>(1 - 0.801888999613917),           // Q = 1 - P
 217       tolerance);
 218    test_spot(
 219       static_cast<RealType>(126),   // degrees of freedom
 220       static_cast<RealType>(-2),   // non centrality
 221       static_cast<RealType>(-4.33),   // T
 222       static_cast<RealType>(1.252846196792878e-2),       // Probability of result (CDF), P
 223       static_cast<RealType>(1 - 1.252846196792878e-2),           // Q = 1 - P
 224       tolerance);
 225    test_spot(
 226       static_cast<RealType>(20),   // degrees of freedom
 227       static_cast<RealType>(23),   // non centrality
 228       static_cast<RealType>(23),   // T
 229       static_cast<RealType>(0.460134400391924),       // Probability of result (CDF), P
 230       static_cast<RealType>(1 - 0.460134400391924),           // Q = 1 - P
 231       tolerance);
 232    test_spot(
 233       static_cast<RealType>(20),   // degrees of freedom
 234       static_cast<RealType>(33),   // non centrality
 235       static_cast<RealType>(34),   // T
 236       static_cast<RealType>(0.532008386378725),       // Probability of result (CDF), P
 237       static_cast<RealType>(1 - 0.532008386378725),           // Q = 1 - P
 238       tolerance);
 239    test_spot(
 240       static_cast<RealType>(12),   // degrees of freedom
 241       static_cast<RealType>(38),   // non centrality
 242       static_cast<RealType>(39),   // T
 243       static_cast<RealType>(0.495868184917805),       // Probability of result (CDF), P
 244       static_cast<RealType>(1 - 0.495868184917805),           // Q = 1 - P
 245       tolerance);
 246    test_spot(
 247       static_cast<RealType>(12),   // degrees of freedom
 248       static_cast<RealType>(39),   // non centrality
 249       static_cast<RealType>(39),   // T
 250       static_cast<RealType>(0.446304024668836),       // Probability of result (CDF), P
 251       static_cast<RealType>(1 - 0.446304024668836),           // Q = 1 - P
 252       tolerance);
 253    test_spot(
 254       static_cast<RealType>(200),   // degrees of freedom
 255       static_cast<RealType>(38),   // non centrality
 256       static_cast<RealType>(39),   // T
 257       static_cast<RealType>(0.666194209961795),       // Probability of result (CDF), P
 258       static_cast<RealType>(1 - 0.666194209961795),           // Q = 1 - P
 259       tolerance);
 260    test_spot(
 261       static_cast<RealType>(200),   // degrees of freedom
 262       static_cast<RealType>(42),   // non centrality
 263       static_cast<RealType>(40),   // T
 264       static_cast<RealType>(0.179292265426085),       // Probability of result (CDF), P
 265       static_cast<RealType>(1 - 0.179292265426085),           // Q = 1 - P
 266       tolerance);
 267
 268    // From https://svn.boost.org/trac/boost/ticket/10480.
 269    // Test value from Mathematica N[CDF[NoncentralStudentTDistribution[2, 4], 5], 35]:
 270    test_spot(
 271       static_cast<RealType>(2),   // degrees of freedom
 272       static_cast<RealType>(4),   // non centrality
 273       static_cast<RealType>(5),   // T
 274       static_cast<RealType>(0.53202069866995310466912357978934321L),       // Probability of result (CDF), P
 275       static_cast<RealType>(1 - 0.53202069866995310466912357978934321L),           // Q = 1 - P
 276       tolerance);
 277
 278    /* This test fails
 279    "Result of tgamma is too large to represent" at naive_mean check for max and infinity.
 280    if (std::numeric_limits<RealType>::has_infinity)
 281    {
 282    test_spot(
 283    //static_cast<RealType>(std::numeric_limits<RealType>::infinity()),   // degrees of freedom
 284    static_cast<RealType>((std::numeric_limits<RealType>::max)()),   // degrees of freedom
 285    static_cast<RealType>(10),   // non centrality
 286    static_cast<RealType>(11),   // T
 287    static_cast<RealType>(0.84134474606854293),       // Probability of result (CDF), P
 288    static_cast<RealType>(0.15865525393145707),           // Q = 1 - P
 289    tolerance);
 290    }
 291    */
 292
 293    boost::math::non_central_t_distribution<RealType> dist(static_cast<RealType>(8), static_cast<RealType>(12));
 294    BOOST_CHECK_CLOSE(pdf(dist, 12), static_cast<RealType>(1.235329715425894935157684607751972713457e-1L), tolerance);
 295    BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution<RealType>(126, -2), -4), static_cast<RealType>(5.797932289365814702402873546466798025787e-2L), tolerance);
 296    BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution<RealType>(126, 2), 4), static_cast<RealType>(5.797932289365814702402873546466798025787e-2L), tolerance);
 297    BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution<RealType>(126, 2), 0), static_cast<RealType>(5.388394890639957139696546086044839573749e-2L), tolerance);
 298
 299    // Error handling checks:
 300    //check_out_of_range<boost::math::non_central_t_distribution<RealType> >(1, 1);  // Fails one check because df for this distribution *can* be infinity.
 301    BOOST_MATH_CHECK_THROW(pdf(boost::math::non_central_t_distribution<RealType>(0, 1), 0), std::domain_error);
 302    BOOST_MATH_CHECK_THROW(pdf(boost::math::non_central_t_distribution<RealType>(-1, 1), 0), std::domain_error);
 303    BOOST_MATH_CHECK_THROW(quantile(boost::math::non_central_t_distribution<RealType>(1, 1), -1), std::domain_error);
 304    BOOST_MATH_CHECK_THROW(quantile(boost::math::non_central_t_distribution<RealType>(1, 1), 2), std::domain_error);
 305 } // template <class RealType>void test_spots(RealType)
 306
 307 template <class T>
 308 T nct_cdf(T df, T nc, T x)
 309 {
 310    return cdf(boost::math::non_central_t_distribution<T>(df, nc), x);
 311 }
 312
 313 template <class T>
 314 T nct_ccdf(T df, T nc, T x)
 315 {
 316    return cdf(complement(boost::math::non_central_t_distribution<T>(df, nc), x));
 317 }
 318
 319 template <typename Real, typename T>
 320 void do_test_nc_t(T& data, const char* type_name, const char* test)
 321 {
 322    typedef Real                   value_type;
 323
 324    std::cout << "Testing: " << test << std::endl;
 325
 326 #ifdef NC_T_CDF_FUNCTION_TO_TEST
 327    value_type(*fp1)(value_type, value_type, value_type) = NC_T_CDF_FUNCTION_TO_TEST;
 328 #else
 329    value_type(*fp1)(value_type, value_type, value_type) = nct_cdf;
 330 #endif
 331    boost::math::tools::test_result<value_type> result;
 332
 333 #if !(defined(ERROR_REPORTING_MODE) && !defined(NC_T_CDF_FUNCTION_TO_TEST))
 334    result = boost::math::tools::test_hetero<Real>(
 335       data,
 336       bind_func<Real>(fp1, 0, 1, 2),
 337       extract_result<Real>(3));
 338    handle_test_result(result, data[result.worst()], result.worst(),
 339       type_name, "non central t CDF", test);
 340 #endif
 341
 342 #if !(defined(ERROR_REPORTING_MODE) && !defined(NC_T_CCDF_FUNCTION_TO_TEST))
 343 #ifdef NC_T_CCDF_FUNCTION_TO_TEST
 344    fp1 = NC_T_CCDF_FUNCTION_TO_TEST;
 345 #else
 346    fp1 = nct_ccdf;
 347 #endif
 348    result = boost::math::tools::test_hetero<Real>(
 349       data,
 350       bind_func<Real>(fp1, 0, 1, 2),
 351       extract_result<Real>(4));
 352    handle_test_result(result, data[result.worst()], result.worst(),
 353       type_name, "non central t CDF complement", test);
 354
 355    std::cout << std::endl;
 356 #endif
 357 }
 358
 359 template <typename Real, typename T>
 360 void quantile_sanity_check(T& data, const char* type_name, const char* test)
 361 {
 362 #ifndef ERROR_REPORTING_MODE
 363    typedef Real                   value_type;
 364
 365    //
 366    // Tests with type real_concept take rather too long to run, so
 367    // for now we'll disable them:
 368    //
 369    if(!boost::is_floating_point<value_type>::value)
 370       return;
 371
 372    std::cout << "Testing: " << type_name << " quantile sanity check, with tests " << test << std::endl;
 373
 374    //
 375    // These sanity checks test for a round trip accuracy of one half
 376    // of the bits in T, unless T is type float, in which case we check
 377    // for just one decimal digit.  The problem here is the sensitivity
 378    // of the functions, not their accuracy.  This test data was generated
 379    // for the forward functions, which means that when it is used as
 380    // the input to the inverses then it is necessarily inexact.  This rounding
 381    // of the input is what makes the data unsuitable for use as an accuracy check,
 382    // and also demonstrates that you can't in general round-trip these functions.
 383    // It is however a useful sanity check.
 384    //
 385    value_type precision = static_cast<value_type>(ldexp(1.0, 1 - boost::math::policies::digits<value_type, boost::math::policies::policy<> >() / 2)) * 100;
 386    if(boost::math::policies::digits<value_type, boost::math::policies::policy<> >() < 50)
 387       precision = 1;   // 1% or two decimal digits, all we can hope for when the input is truncated to float
 388
 389    for(unsigned i = 0; i < data.size(); ++i)
 390    {
 391       if(data[i][3] == 0)
 392       {
 393          BOOST_CHECK(0 == quantile(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), data[i][3]));
 394       }
 395       else if(data[i][3] < 0.9999f)
 396       {
 397          value_type p = quantile(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), data[i][3]);
 398          value_type pt = data[i][2];
 399          BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
 400       }
 401       if(data[i][4] == 0)
 402       {
 403          BOOST_CHECK(0 == quantile(complement(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), data[i][3])));
 404       }
 405       else if(data[i][4] < 0.9999f)
 406       {
 407          value_type p = quantile(complement(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), data[i][4]));
 408          value_type pt = data[i][2];
 409          BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
 410       }
 411       if(boost::math::tools::digits<value_type>() > 50)
 412       {
 413          //
 414          // Sanity check mode, the accuracy of
 415          // the mode is at *best* the square root of the accuracy of the PDF:
 416          //
 417 #ifndef BOOST_NO_EXCEPTIONS
 418          try{
 419             value_type m = mode(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]));
 420             value_type p = pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m);
 421             value_type delta = (std::max)(fabs(m * sqrt(precision) * 50), sqrt(precision) * 50);
 422             BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m + delta) <= p, i);
 423             BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m - delta) <= p, i);
 424          }
 425          catch(const boost::math::evaluation_error&) {}
 426 #endif
 427 #if 0
 428          //
 429          // Sanity check degrees-of-freedom finder, don't bother at float
 430          // precision though as there's not enough data in the probability
 431          // values to get back to the correct degrees of freedom or
 432          // non-centrality parameter:
 433          //
 434          try{
 435             if((data[i][3] < 0.99) && (data[i][3] != 0))
 436             {
 437                BOOST_CHECK_CLOSE_EX(
 438                   boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], data[i][3]),
 439                   data[i][0], precision, i);
 440                BOOST_CHECK_CLOSE_EX(
 441                   boost::math::non_central_t_distribution<value_type>::find_non_centrality(data[i][0], data[i][2], data[i][3]),
 442                   data[i][1], precision, i);
 443             }
 444             if((data[i][4] < 0.99) && (data[i][4] != 0))
 445             {
 446                BOOST_CHECK_CLOSE_EX(
 447                   boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], data[i][4])),
 448                   data[i][0], precision, i);
 449                BOOST_CHECK_CLOSE_EX(
 450                   boost::math::non_central_t_distribution<value_type>::find_non_centrality(boost::math::complement(data[i][0], data[i][2], data[i][4])),
 451                   data[i][1], precision, i);
 452             }
 453          }
 454          catch(const std::exception& e)
 455          {
 456             BOOST_ERROR(e.what());
 457          }
 458 #endif
 459       }
 460    }
 461 #endif
 462 }
 463
 464 template <typename T>
 465 void test_accuracy(T, const char* type_name)
 466 {
 467 #include "nct.ipp"
 468    do_test_nc_t<T>(nct, type_name, "Non Central T");
 469    quantile_sanity_check<T>(nct, type_name, "Non Central T");
 470    if(std::numeric_limits<T>::is_specialized)
 471    {
 472       //
 473       // Don't run these tests for real_concept: they take too long and don't converge
 474       // without numeric_limits and lanczos support:
 475       //
 476 #include "nct_small_delta.ipp"
 477       do_test_nc_t<T>(nct_small_delta, type_name, "Non Central T (small non-centrality)");
 478       quantile_sanity_check<T>(nct_small_delta, type_name, "Non Central T (small non-centrality)");
 479 #include "nct_asym.ipp"
 480       do_test_nc_t<T>(nct_asym, type_name, "Non Central T (large parameters)");
 481       quantile_sanity_check<T>(nct_asym, type_name, "Non Central T (large parameters)");
 482    }
 483 }
 484
 485
 486 template <class RealType>
 487 void test_big_df(RealType)
 488 {
 489    using namespace boost::math;
 490
 491    if(typeid(RealType) != typeid(boost::math::concepts::real_concept))
 492    { // Ordinary floats only.
 493       // Could also test if (std::numeric_limits<RealType>::is_specialized);
 494
 495       RealType tolerance = 10 * boost::math::tools::epsilon<RealType>(); // static_cast<RealType>(1e-14); //
 496       std::cout.precision(17); // Note: need to reset after calling BOOST_CHECK_s
 497       // due to buglet in Boost.test that fails to restore precision corrrectly.
 498
 499       // Test for large degrees of freedom when should be same as normal.
 500       RealType inf =
 501          (std::numeric_limits<RealType>::has_infinity) ?
 502          std::numeric_limits<RealType>::infinity()
 503          :
 504          boost::math::tools::max_value<RealType>();
 505       RealType nan = std::numeric_limits<RealType>::quiet_NaN();
 506
 507       // Tests for df = max_value and infinity.
 508       RealType max_val = boost::math::tools::max_value<RealType>();
 509       non_central_t_distribution<RealType> maxdf(max_val, 0);
 510       BOOST_CHECK_EQUAL(maxdf.degrees_of_freedom(), max_val);
 511
 512       non_central_t_distribution<RealType> infdf(inf, 0);
 513       BOOST_CHECK_EQUAL(infdf.degrees_of_freedom(), inf);
 514       BOOST_CHECK_EQUAL(mean(infdf), 0);
 515       BOOST_CHECK_EQUAL(mean(maxdf), 0);
 516       BOOST_CHECK_EQUAL(variance(infdf), 1);
 517       BOOST_CHECK_EQUAL(variance(maxdf), 1);
 518       BOOST_CHECK_EQUAL(skewness(infdf), 0);
 519       BOOST_CHECK_EQUAL(skewness(maxdf), 0);
 520       BOOST_CHECK_EQUAL(kurtosis_excess(infdf), 3);
 521       BOOST_CHECK_CLOSE_FRACTION(kurtosis_excess(maxdf), static_cast<RealType>(3), tolerance);
 522
 523       // Bad df examples.
 524 #ifndef BOOST_NO_EXCEPTIONS
 525       BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType> minfdf(-inf, 0), std::domain_error);
 526       BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType> minfdf(nan, 0), std::domain_error);
 527       BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType> minfdf(-nan, 0), std::domain_error);
 528 #else
 529       BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType>(-inf, 0), std::domain_error);
 530       BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType>(nan, 0), std::domain_error);
 531       BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType>(-nan, 0), std::domain_error);
 532 #endif
 533
 534
 535       // BOOST_CHECK_CLOSE_FRACTION(pdf(infdf, 0), static_cast<RealType>(0.3989422804014326779399460599343818684759L), tolerance);
 536       BOOST_CHECK_CLOSE_FRACTION(pdf(maxdf, 0), boost::math::constants::one_div_root_two_pi<RealType>(), tolerance);
 537       BOOST_CHECK_CLOSE_FRACTION(pdf(infdf, 0), boost::math::constants::one_div_root_two_pi<RealType>(), tolerance);
 538       BOOST_CHECK_CLOSE_FRACTION(cdf(infdf, 0), boost::math::constants::half<RealType>(), tolerance);
 539       BOOST_CHECK_CLOSE_FRACTION(cdf(maxdf, 0), boost::math::constants::half<RealType>(), tolerance);
 540
 541       // non-centrality delta = 10
 542       // Degrees of freedom = Max value and  = infinity should be very close.
 543       non_central_t_distribution<RealType> maxdf10(max_val, 10);
 544       non_central_t_distribution<RealType> infdf10(inf, 10);
 545       BOOST_CHECK_EQUAL(infdf10.degrees_of_freedom(), inf);
 546       BOOST_CHECK_EQUAL(infdf10.non_centrality(), 10);
 547       BOOST_CHECK_EQUAL(mean(infdf10), 10);
 548       BOOST_CHECK_CLOSE_FRACTION(mean(maxdf10), static_cast<RealType>(10), tolerance);
 549
 550       BOOST_CHECK_CLOSE_FRACTION(pdf(infdf10, 11), pdf(maxdf10, 11), tolerance); //
 551
 552       BOOST_CHECK_CLOSE_FRACTION(cdf(complement(infdf10, 11)), 1 - cdf(infdf10, 11), tolerance); //
 553       BOOST_CHECK_CLOSE_FRACTION(cdf(complement(maxdf10, 11)), 1 - cdf(maxdf10, 11), tolerance); //
 554       BOOST_CHECK_CLOSE_FRACTION(cdf(complement(infdf10, 11)), 1 - cdf(maxdf10, 11), tolerance); //
 555       std::cout.precision(17);
 556       //std::cout  << "cdf(maxdf10, 11)  = " << cdf(maxdf10, 11) << ' ' << cdf(complement(maxdf10, 11)) << endl;
 557       //std::cout  << "cdf(infdf10, 11)  = " << cdf(infdf10, 11) << ' ' << cdf(complement(infdf10, 11)) << endl;
 558       //std::cout  << "quantile(maxdf10, 0.5)  = " << quantile(maxdf10, 0.5) << std::endl; // quantile(maxdf10, 0.5)  = 10.000000000000004
 559       //std::cout  << "quantile(infdf10, 0.5) = " << ' ' << quantile(infdf10, 0.5) << std::endl; // quantile(infdf10, 0.5) =  10
 560
 561       BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.5), static_cast<RealType>(10), tolerance);
 562       BOOST_CHECK_CLOSE_FRACTION(quantile(maxdf10, 0.5), static_cast<RealType>(10), tolerance);
 563
 564       BOOST_TEST_MESSAGE("non_central_t_distribution<RealType> infdf100(inf, 100);");
 565       non_central_t_distribution<RealType> infdf100(inf, 100);
 566       BOOST_TEST_MESSAGE("non_central_t_distribution<RealType> maxdf100(max_val, 100);");
 567       non_central_t_distribution<RealType> maxdf100(max_val, 100);
 568       BOOST_TEST_MESSAGE("BOOST_CHECK_CLOSE_FRACTION(quantile(infdf100, 0.5), static_cast<RealType>(100), tolerance);");
 569       BOOST_CHECK_CLOSE_FRACTION(quantile(infdf100, 0.5), static_cast<RealType>(100), tolerance);
 570       BOOST_TEST_MESSAGE("BOOST_CHECK_CLOSE_FRACTION(quantile(maxdf100, 0.5), static_cast<RealType>(100), tolerance);");
 571       BOOST_CHECK_CLOSE_FRACTION(quantile(maxdf100, 0.5), static_cast<RealType>(100), tolerance);
 572       { // Loop back.
 573          RealType p = static_cast<RealType>(0.01);
 574          RealType x = quantile(infdf10, p);
 575          RealType c = cdf(infdf10, x);
 576          BOOST_CHECK_CLOSE_FRACTION(c, p, tolerance);
 577       }
 578     {
 579        RealType q = static_cast<RealType>(0.99);
 580        RealType x = quantile(complement(infdf10, q));
 581        RealType c = cdf(complement(infdf10, x));
 582        BOOST_CHECK_CLOSE_FRACTION(c, q, tolerance);
 583     }
 584     { // Loop back.
 585        RealType p = static_cast<RealType>(0.99);
 586        RealType x = quantile(infdf10, p);
 587        RealType c = cdf(infdf10, x);
 588        BOOST_CHECK_CLOSE_FRACTION(c, p, tolerance);
 589     }
 590     {
 591        RealType q = static_cast<RealType>(0.01);
 592        RealType x = quantile(complement(infdf10, q));
 593        RealType c = cdf(complement(infdf10, x));
 594        BOOST_CHECK_CLOSE_FRACTION(c, q, tolerance * 2); // c{0.0100000128} and q{0.00999999978}
 595     }
 596
 597     //RealType cinf = quantile(infdf10, 0.25);
 598     //std::cout << cinf << ' ' << cdf(infdf10, cinf) << std::endl; // 9.32551 0.25
 599
 600     //RealType cmax = quantile(maxdf10, 0.25);
 601     //std::cout << cmax << ' ' << cdf(maxdf10, cmax) << std::endl; //  9.32551 0.25
 602
 603     //RealType cinfc = quantile(complement(infdf10, 0.75));
 604     //std::cout << cinfc << ' ' << cdf(infdf10, cinfc) << std::endl; // 9.32551 0.25
 605
 606     //RealType cmaxc = quantile(complement(maxdf10, 0.75));
 607     //std::cout << cmaxc << ' ' << cdf(maxdf10, cmaxc) << std::endl; // 9.32551 0.25
 608
 609     BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.5), quantile(maxdf10, 0.5), tolerance); //
 610     BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.2), quantile(maxdf10, 0.2), tolerance); //
 611     BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.8), quantile(maxdf10, 0.8), tolerance); //
 612
 613     BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.25), quantile(complement(infdf10, 0.75)), tolerance); //
 614     BOOST_CHECK_CLOSE_FRACTION(quantile(complement(infdf10, 0.5)), quantile(complement(maxdf10, 0.5)), tolerance); //
 615
 616     BOOST_CHECK_CLOSE_FRACTION(quantile(maxdf10, 0.25), quantile(complement(maxdf10, 0.75)), tolerance); //
 617
 618     BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.99), quantile(complement(infdf10, 0.01)), tolerance); //
 619     BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.4), quantile(complement(infdf10, 0.6)), tolerance); //
 620     BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.01), quantile(complement(infdf10, 1 - 0.01)), tolerance); //
 621    }
 622 } // void test_big_df(RealType)
 623
 624 template <class RealType>
 625 void test_ignore_policy(RealType)
 626 {
 627    // Check on returns when errors are ignored.
 628    if((typeid(RealType) != typeid(boost::math::concepts::real_concept))
 629       && std::numeric_limits<RealType>::has_infinity
 630       && std::numeric_limits<RealType>::has_quiet_NaN
 631       )
 632    { // Ordinary floats only.
 633
 634       using namespace boost::math;
 635       //   RealType inf = std::numeric_limits<RealType>::infinity();
 636       RealType nan = std::numeric_limits<RealType>::quiet_NaN();
 637
 638       using boost::math::policies::policy;
 639       // Types of error whose action can be altered by policies:.
 640       //using boost::math::policies::evaluation_error;
 641       //using boost::math::policies::domain_error;
 642       //using boost::math::policies::overflow_error;
 643       //using boost::math::policies::underflow_error;
 644       //using boost::math::policies::domain_error;
 645       //using boost::math::policies::pole_error;
 646
 647       //// Actions on error (in enum error_policy_type):
 648       //using boost::math::policies::errno_on_error;
 649       //using boost::math::policies::ignore_error;
 650       //using boost::math::policies::throw_on_error;
 651       //using boost::math::policies::denorm_error;
 652       //using boost::math::policies::pole_error;
 653       //using boost::math::policies::user_error;
 654
 655       typedef policy<
 656          boost::math::policies::domain_error<boost::math::policies::ignore_error>,
 657          boost::math::policies::overflow_error<boost::math::policies::ignore_error>,
 658          boost::math::policies::underflow_error<boost::math::policies::ignore_error>,
 659          boost::math::policies::denorm_error<boost::math::policies::ignore_error>,
 660          boost::math::policies::pole_error<boost::math::policies::ignore_error>,
 661          boost::math::policies::evaluation_error<boost::math::policies::ignore_error>
 662       > ignore_all_policy;
 663
 664       typedef non_central_t_distribution<RealType, ignore_all_policy> ignore_error_non_central_t;
 665
 666       // Only test NaN and infinity if type has these features (realconcept returns zero).
 667       // Integers are always converted to RealType,
 668       // others requires static cast to RealType from long double.
 669
 670       if(std::numeric_limits<RealType>::has_quiet_NaN)
 671       {
 672          // Mean
 673          BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(-nan, 0))));
 674          BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(+nan, 0))));
 675          BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(-1, 0))));
 676          BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(0, 0))));
 677          BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(1, 0))));
 678          BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(2, nan))));
 679          BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(nan, nan))));
 680          BOOST_CHECK(boost::math::isfinite(mean(ignore_error_non_central_t(2, 0)))); // OK
 681
 682          // Variance
 683          BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(nan, 0))));
 684          BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(1, nan))));
 685          BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(2, nan))));
 686          BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(-1, 0))));
 687          BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(0, 0))));
 688          BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(1, 0))));
 689          BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(static_cast<RealType>(1.7L), 0))));
 690          BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(2, 0))));
 691
 692          // Skewness
 693          BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(std::numeric_limits<RealType>::quiet_NaN(), 0))));
 694          BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(-1, 0))));
 695          BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(0, 0))));
 696          BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(1, 0))));
 697          BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(2, 0))));
 698          BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(3, 0))));
 699
 700          // Kurtosis
 701          BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(std::numeric_limits<RealType>::quiet_NaN(), 0))));
 702          BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(-1, 0))));
 703          BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(0, 0))));
 704          BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(1, 0))));
 705          BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(2, 0))));
 706          BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(static_cast<RealType>(2.0001L), 0))));
 707          BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(3, 0))));
 708          BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(4, 0))));
 709
 710          // Kurtosis excess
 711          BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(std::numeric_limits<RealType>::quiet_NaN(), 0))));
 712          BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(-1, 0))));
 713          BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(0, 0))));
 714          BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(1, 0))));
 715          BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(2, 0))));
 716          BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(static_cast<RealType>(2.0001L), 0))));
 717          BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(3, 0))));
 718          BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(4, 0))));
 719       } // has_quiet_NaN
 720       BOOST_CHECK(boost::math::isfinite(mean(ignore_error_non_central_t(1 + std::numeric_limits<RealType>::epsilon(), 0))));
 721       BOOST_CHECK(boost::math::isfinite(variance(ignore_error_non_central_t(2 + 2 * std::numeric_limits<RealType>::epsilon(), 0))));
 722       BOOST_CHECK(boost::math::isfinite(variance(ignore_error_non_central_t(static_cast<RealType>(2.0001L), 0))));
 723       BOOST_CHECK(boost::math::isfinite(variance(ignore_error_non_central_t(2 + 2 * std::numeric_limits<RealType>::epsilon(), 0))));
 724       BOOST_CHECK(boost::math::isfinite(skewness(ignore_error_non_central_t(3 + 3 * std::numeric_limits<RealType>::epsilon(), 0))));
 725       BOOST_CHECK(boost::math::isfinite(kurtosis(ignore_error_non_central_t(4 + 4 * std::numeric_limits<RealType>::epsilon(), 0))));
 726       BOOST_CHECK(boost::math::isfinite(kurtosis(ignore_error_non_central_t(static_cast<RealType>(4.0001L), 0))));
 727
 728       // check_out_of_range<non_central_t_distribution<RealType> >(1, 0); // Fails one check because allows df = infinity.
 729       check_support<non_central_t_distribution<RealType> >(non_central_t_distribution<RealType>(1, 0));
 730    } // ordinary floats.
 731 } // template <class RealType> void test_ignore_policy(RealType)