ceph/src/boost/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/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 typename T::value_type row_type;
 323    typedef Real                   value_type;
 324
 325    std::cout << "Testing: " << test << std::endl;
 326
 327 #ifdef NC_T_CDF_FUNCTION_TO_TEST
 328    value_type(*fp1)(value_type, value_type, value_type) = NC_T_CDF_FUNCTION_TO_TEST;
 329 #else
 330    value_type(*fp1)(value_type, value_type, value_type) = nct_cdf;
 331 #endif
 332    boost::math::tools::test_result<value_type> result;
 333
 334 #if !(defined(ERROR_REPORTING_MODE) && !defined(NC_T_CDF_FUNCTION_TO_TEST))
 335    result = boost::math::tools::test_hetero<Real>(
 336       data,
 337       bind_func<Real>(fp1, 0, 1, 2),
 338       extract_result<Real>(3));
 339    handle_test_result(result, data[result.worst()], result.worst(),
 340       type_name, "non central t CDF", test);
 341 #endif
 342
 343 #if !(defined(ERROR_REPORTING_MODE) && !defined(NC_T_CCDF_FUNCTION_TO_TEST))
 344 #ifdef NC_T_CCDF_FUNCTION_TO_TEST
 345    fp1 = NC_T_CCDF_FUNCTION_TO_TEST;
 346 #else
 347    fp1 = nct_ccdf;
 348 #endif
 349    result = boost::math::tools::test_hetero<Real>(
 350       data,
 351       bind_func<Real>(fp1, 0, 1, 2),
 352       extract_result<Real>(4));
 353    handle_test_result(result, data[result.worst()], result.worst(),
 354       type_name, "non central t CDF complement", test);
 355
 356    std::cout << std::endl;
 357 #endif
 358 }
 359
 360 template <typename Real, typename T>
 361 void quantile_sanity_check(T& data, const char* type_name, const char* test)
 362 {
 363 #ifndef ERROR_REPORTING_MODE
 364    typedef typename T::value_type row_type;
 365    typedef Real                   value_type;
 366
 367    //
 368    // Tests with type real_concept take rather too long to run, so
 369    // for now we'll disable them:
 370    //
 371    if(!boost::is_floating_point<value_type>::value)
 372       return;
 373
 374    std::cout << "Testing: " << type_name << " quantile sanity check, with tests " << test << std::endl;
 375
 376    //
 377    // These sanity checks test for a round trip accuracy of one half
 378    // of the bits in T, unless T is type float, in which case we check
 379    // for just one decimal digit.  The problem here is the sensitivity
 380    // of the functions, not their accuracy.  This test data was generated
 381    // for the forward functions, which means that when it is used as
 382    // the input to the inverses then it is necessarily inexact.  This rounding
 383    // of the input is what makes the data unsuitable for use as an accuracy check,
 384    // and also demonstrates that you can't in general round-trip these functions.
 385    // It is however a useful sanity check.
 386    //
 387    value_type precision = static_cast<value_type>(ldexp(1.0, 1 - boost::math::policies::digits<value_type, boost::math::policies::policy<> >() / 2)) * 100;
 388    if(boost::math::policies::digits<value_type, boost::math::policies::policy<> >() < 50)
 389       precision = 1;   // 1% or two decimal digits, all we can hope for when the input is truncated to float
 390
 391    for(unsigned i = 0; i < data.size(); ++i)
 392    {
 393       if(data[i][3] == 0)
 394       {
 395          BOOST_CHECK(0 == quantile(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), data[i][3]));
 396       }
 397       else if(data[i][3] < 0.9999f)
 398       {
 399          value_type p = quantile(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), data[i][3]);
 400          value_type pt = data[i][2];
 401          BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
 402       }
 403       if(data[i][4] == 0)
 404       {
 405          BOOST_CHECK(0 == quantile(complement(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), data[i][3])));
 406       }
 407       else if(data[i][4] < 0.9999f)
 408       {
 409          value_type p = quantile(complement(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), data[i][4]));
 410          value_type pt = data[i][2];
 411          BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
 412       }
 413       if(boost::math::tools::digits<value_type>() > 50)
 414       {
 415          //
 416          // Sanity check mode, the accuracy of
 417          // the mode is at *best* the square root of the accuracy of the PDF:
 418          //
 419 #ifndef BOOST_NO_EXCEPTIONS
 420          try{
 421             value_type m = mode(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]));
 422             value_type p = pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m);
 423             value_type delta = (std::max)(fabs(m * sqrt(precision) * 50), sqrt(precision) * 50);
 424             BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m + delta) <= p, i);
 425             BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m - delta) <= p, i);
 426          }
 427          catch(const boost::math::evaluation_error&) {}
 428 #endif
 429 #if 0
 430          //
 431          // Sanity check degrees-of-freedom finder, don't bother at float
 432          // precision though as there's not enough data in the probability
 433          // values to get back to the correct degrees of freedom or
 434          // non-centrality parameter:
 435          //
 436          try{
 437             if((data[i][3] < 0.99) && (data[i][3] != 0))
 438             {
 439                BOOST_CHECK_CLOSE_EX(
 440                   boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], data[i][3]),
 441                   data[i][0], precision, i);
 442                BOOST_CHECK_CLOSE_EX(
 443                   boost::math::non_central_t_distribution<value_type>::find_non_centrality(data[i][0], data[i][2], data[i][3]),
 444                   data[i][1], precision, i);
 445             }
 446             if((data[i][4] < 0.99) && (data[i][4] != 0))
 447             {
 448                BOOST_CHECK_CLOSE_EX(
 449                   boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], data[i][4])),
 450                   data[i][0], precision, i);
 451                BOOST_CHECK_CLOSE_EX(
 452                   boost::math::non_central_t_distribution<value_type>::find_non_centrality(boost::math::complement(data[i][0], data[i][2], data[i][4])),
 453                   data[i][1], precision, i);
 454             }
 455          }
 456          catch(const std::exception& e)
 457          {
 458             BOOST_ERROR(e.what());
 459          }
 460 #endif
 461       }
 462    }
 463 #endif
 464 }
 465
 466 template <typename T>
 467 void test_accuracy(T, const char* type_name)
 468 {
 469 #include "nct.ipp"
 470    do_test_nc_t<T>(nct, type_name, "Non Central T");
 471    quantile_sanity_check<T>(nct, type_name, "Non Central T");
 472    if(std::numeric_limits<T>::is_specialized)
 473    {
 474       //
 475       // Don't run these tests for real_concept: they take too long and don't converge
 476       // without numeric_limits and lanczos support:
 477       //
 478 #include "nct_small_delta.ipp"
 479       do_test_nc_t<T>(nct_small_delta, type_name, "Non Central T (small non-centrality)");
 480       quantile_sanity_check<T>(nct_small_delta, type_name, "Non Central T (small non-centrality)");
 481 #include "nct_asym.ipp"
 482       do_test_nc_t<T>(nct_asym, type_name, "Non Central T (large parameters)");
 483       quantile_sanity_check<T>(nct_asym, type_name, "Non Central T (large parameters)");
 484    }
 485 }
 486
 487
 488 template <class RealType>
 489 void test_big_df(RealType)
 490 {
 491    using namespace boost::math;
 492
 493    if(typeid(RealType) != typeid(boost::math::concepts::real_concept))
 494    { // Ordinary floats only.
 495       // Could also test if (std::numeric_limits<RealType>::is_specialized);
 496
 497       RealType tolerance = 10 * boost::math::tools::epsilon<RealType>(); // static_cast<RealType>(1e-14); //
 498       std::cout.precision(17); // Note: need to reset after calling BOOST_CHECK_s
 499       // due to buglet in Boost.test that fails to restore precision corrrectly.
 500
 501       // Test for large degrees of freedom when should be same as normal.
 502       RealType inf =
 503          (std::numeric_limits<RealType>::has_infinity) ?
 504          std::numeric_limits<RealType>::infinity()
 505          :
 506          boost::math::tools::max_value<RealType>();
 507       RealType nan = std::numeric_limits<RealType>::quiet_NaN();
 508
 509       // Tests for df = max_value and infinity.
 510       RealType max_val = boost::math::tools::max_value<RealType>();
 511       non_central_t_distribution<RealType> maxdf(max_val, 0);
 512       BOOST_CHECK_EQUAL(maxdf.degrees_of_freedom(), max_val);
 513
 514       non_central_t_distribution<RealType> infdf(inf, 0);
 515       BOOST_CHECK_EQUAL(infdf.degrees_of_freedom(), inf);
 516       BOOST_CHECK_EQUAL(mean(infdf), 0);
 517       BOOST_CHECK_EQUAL(mean(maxdf), 0);
 518       BOOST_CHECK_EQUAL(variance(infdf), 1);
 519       BOOST_CHECK_EQUAL(variance(maxdf), 1);
 520       BOOST_CHECK_EQUAL(skewness(infdf), 0);
 521       BOOST_CHECK_EQUAL(skewness(maxdf), 0);
 522       BOOST_CHECK_EQUAL(kurtosis_excess(infdf), 3);
 523       BOOST_CHECK_CLOSE_FRACTION(kurtosis_excess(maxdf), static_cast<RealType>(3), tolerance);
 524
 525       // Bad df examples.
 526 #ifndef BOOST_NO_EXCEPTIONS
 527       BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType> minfdf(-inf, 0), std::domain_error);
 528       BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType> minfdf(nan, 0), std::domain_error);
 529       BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType> minfdf(-nan, 0), std::domain_error);
 530 #else
 531       BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType>(-inf, 0), std::domain_error);
 532       BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType>(nan, 0), std::domain_error);
 533       BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType>(-nan, 0), std::domain_error);
 534 #endif
 535
 536
 537       // BOOST_CHECK_CLOSE_FRACTION(pdf(infdf, 0), static_cast<RealType>(0.3989422804014326779399460599343818684759L), tolerance);
 538       BOOST_CHECK_CLOSE_FRACTION(pdf(maxdf, 0), boost::math::constants::one_div_root_two_pi<RealType>(), tolerance);
 539       BOOST_CHECK_CLOSE_FRACTION(pdf(infdf, 0), boost::math::constants::one_div_root_two_pi<RealType>(), tolerance);
 540       BOOST_CHECK_CLOSE_FRACTION(cdf(infdf, 0), boost::math::constants::half<RealType>(), tolerance);
 541       BOOST_CHECK_CLOSE_FRACTION(cdf(maxdf, 0), boost::math::constants::half<RealType>(), tolerance);
 542
 543       // non-centrality delta = 10
 544       // Degrees of freedom = Max value and  = infinity should be very close.
 545       non_central_t_distribution<RealType> maxdf10(max_val, 10);
 546       non_central_t_distribution<RealType> infdf10(inf, 10);
 547       BOOST_CHECK_EQUAL(infdf10.degrees_of_freedom(), inf);
 548       BOOST_CHECK_EQUAL(infdf10.non_centrality(), 10);
 549       BOOST_CHECK_EQUAL(mean(infdf10), 10);
 550       BOOST_CHECK_CLOSE_FRACTION(mean(maxdf10), static_cast<RealType>(10), tolerance);
 551
 552       BOOST_CHECK_CLOSE_FRACTION(pdf(infdf10, 11), pdf(maxdf10, 11), tolerance); //
 553
 554       BOOST_CHECK_CLOSE_FRACTION(cdf(complement(infdf10, 11)), 1 - cdf(infdf10, 11), tolerance); //
 555       BOOST_CHECK_CLOSE_FRACTION(cdf(complement(maxdf10, 11)), 1 - cdf(maxdf10, 11), tolerance); //
 556       BOOST_CHECK_CLOSE_FRACTION(cdf(complement(infdf10, 11)), 1 - cdf(maxdf10, 11), tolerance); //
 557       std::cout.precision(17);
 558       //std::cout  << "cdf(maxdf10, 11)  = " << cdf(maxdf10, 11) << ' ' << cdf(complement(maxdf10, 11)) << endl;
 559       //std::cout  << "cdf(infdf10, 11)  = " << cdf(infdf10, 11) << ' ' << cdf(complement(infdf10, 11)) << endl;
 560       //std::cout  << "quantile(maxdf10, 0.5)  = " << quantile(maxdf10, 0.5) << std::endl; // quantile(maxdf10, 0.5)  = 10.000000000000004
 561       //std::cout  << "quantile(infdf10, 0.5) = " << ' ' << quantile(infdf10, 0.5) << std::endl; // quantile(infdf10, 0.5) =  10
 562
 563       BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.5), static_cast<RealType>(10), tolerance);
 564       BOOST_CHECK_CLOSE_FRACTION(quantile(maxdf10, 0.5), static_cast<RealType>(10), tolerance);
 565
 566       BOOST_TEST_MESSAGE("non_central_t_distribution<RealType> infdf100(inf, 100);");
 567       non_central_t_distribution<RealType> infdf100(inf, 100);
 568       BOOST_TEST_MESSAGE("non_central_t_distribution<RealType> maxdf100(max_val, 100);");
 569       non_central_t_distribution<RealType> maxdf100(max_val, 100);
 570       BOOST_TEST_MESSAGE("BOOST_CHECK_CLOSE_FRACTION(quantile(infdf100, 0.5), static_cast<RealType>(100), tolerance);");
 571       BOOST_CHECK_CLOSE_FRACTION(quantile(infdf100, 0.5), static_cast<RealType>(100), tolerance);
 572       BOOST_TEST_MESSAGE("BOOST_CHECK_CLOSE_FRACTION(quantile(maxdf100, 0.5), static_cast<RealType>(100), tolerance);");
 573       BOOST_CHECK_CLOSE_FRACTION(quantile(maxdf100, 0.5), static_cast<RealType>(100), tolerance);
 574       { // Loop back.
 575          RealType p = static_cast<RealType>(0.01);
 576          RealType x = quantile(infdf10, p);
 577          RealType c = cdf(infdf10, x);
 578          BOOST_CHECK_CLOSE_FRACTION(c, p, tolerance);
 579       }
 580     {
 581        RealType q = static_cast<RealType>(0.99);
 582        RealType x = quantile(complement(infdf10, q));
 583        RealType c = cdf(complement(infdf10, x));
 584        BOOST_CHECK_CLOSE_FRACTION(c, q, tolerance);
 585     }
 586     { // Loop back.
 587        RealType p = static_cast<RealType>(0.99);
 588        RealType x = quantile(infdf10, p);
 589        RealType c = cdf(infdf10, x);
 590        BOOST_CHECK_CLOSE_FRACTION(c, p, tolerance);
 591     }
 592     {
 593        RealType q = static_cast<RealType>(0.01);
 594        RealType x = quantile(complement(infdf10, q));
 595        RealType c = cdf(complement(infdf10, x));
 596        BOOST_CHECK_CLOSE_FRACTION(c, q, tolerance * 2); // c{0.0100000128} and q{0.00999999978}
 597     }
 598
 599     //RealType cinf = quantile(infdf10, 0.25);
 600     //std::cout << cinf << ' ' << cdf(infdf10, cinf) << std::endl; // 9.32551 0.25
 601
 602     //RealType cmax = quantile(maxdf10, 0.25);
 603     //std::cout << cmax << ' ' << cdf(maxdf10, cmax) << std::endl; //  9.32551 0.25
 604
 605     //RealType cinfc = quantile(complement(infdf10, 0.75));
 606     //std::cout << cinfc << ' ' << cdf(infdf10, cinfc) << std::endl; // 9.32551 0.25
 607
 608     //RealType cmaxc = quantile(complement(maxdf10, 0.75));
 609     //std::cout << cmaxc << ' ' << cdf(maxdf10, cmaxc) << std::endl; // 9.32551 0.25
 610
 611     BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.5), quantile(maxdf10, 0.5), tolerance); //
 612     BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.2), quantile(maxdf10, 0.2), tolerance); //
 613     BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.8), quantile(maxdf10, 0.8), tolerance); //
 614
 615     BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.25), quantile(complement(infdf10, 0.75)), tolerance); //
 616     BOOST_CHECK_CLOSE_FRACTION(quantile(complement(infdf10, 0.5)), quantile(complement(maxdf10, 0.5)), tolerance); //
 617
 618     BOOST_CHECK_CLOSE_FRACTION(quantile(maxdf10, 0.25), quantile(complement(maxdf10, 0.75)), tolerance); //
 619
 620     BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.99), quantile(complement(infdf10, 0.01)), tolerance); //
 621     BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.4), quantile(complement(infdf10, 0.6)), tolerance); //
 622     BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.01), quantile(complement(infdf10, 1 - 0.01)), tolerance); //
 623    }
 624 } // void test_big_df(RealType)
 625
 626 template <class RealType>
 627 void test_ignore_policy(RealType)
 628 {
 629    // Check on returns when errors are ignored.
 630    if((typeid(RealType) != typeid(boost::math::concepts::real_concept))
 631       && std::numeric_limits<RealType>::has_infinity
 632       && std::numeric_limits<RealType>::has_quiet_NaN
 633       )
 634    { // Ordinary floats only.
 635
 636       using namespace boost::math;
 637       //   RealType inf = std::numeric_limits<RealType>::infinity();
 638       RealType nan = std::numeric_limits<RealType>::quiet_NaN();
 639
 640       using boost::math::policies::policy;
 641       // Types of error whose action can be altered by policies:.
 642       //using boost::math::policies::evaluation_error;
 643       //using boost::math::policies::domain_error;
 644       //using boost::math::policies::overflow_error;
 645       //using boost::math::policies::underflow_error;
 646       //using boost::math::policies::domain_error;
 647       //using boost::math::policies::pole_error;
 648
 649       //// Actions on error (in enum error_policy_type):
 650       //using boost::math::policies::errno_on_error;
 651       //using boost::math::policies::ignore_error;
 652       //using boost::math::policies::throw_on_error;
 653       //using boost::math::policies::denorm_error;
 654       //using boost::math::policies::pole_error;
 655       //using boost::math::policies::user_error;
 656
 657       typedef policy<
 658          boost::math::policies::domain_error<boost::math::policies::ignore_error>,
 659          boost::math::policies::overflow_error<boost::math::policies::ignore_error>,
 660          boost::math::policies::underflow_error<boost::math::policies::ignore_error>,
 661          boost::math::policies::denorm_error<boost::math::policies::ignore_error>,
 662          boost::math::policies::pole_error<boost::math::policies::ignore_error>,
 663          boost::math::policies::evaluation_error<boost::math::policies::ignore_error>
 664       > ignore_all_policy;
 665
 666       typedef non_central_t_distribution<RealType, ignore_all_policy> ignore_error_non_central_t;
 667
 668       // Only test NaN and infinity if type has these features (realconcept returns zero).
 669       // Integers are always converted to RealType,
 670       // others requires static cast to RealType from long double.
 671
 672       if(std::numeric_limits<RealType>::has_quiet_NaN)
 673       {
 674          // Mean
 675          BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(-nan, 0))));
 676          BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(+nan, 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(0, 0))));
 679          BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(1, 0))));
 680          BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(2, nan))));
 681          BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(nan, nan))));
 682          BOOST_CHECK(boost::math::isfinite(mean(ignore_error_non_central_t(2, 0)))); // OK
 683
 684          // Variance
 685          BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(nan, 0))));
 686          BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(1, nan))));
 687          BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(2, nan))));
 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(0, 0))));
 690          BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(1, 0))));
 691          BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(static_cast<RealType>(1.7L), 0))));
 692          BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(2, 0))));
 693
 694          // Skewness
 695          BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(std::numeric_limits<RealType>::quiet_NaN(), 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(0, 0))));
 698          BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(1, 0))));
 699          BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(2, 0))));
 700          BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(3, 0))));
 701
 702          // Kurtosis
 703          BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(std::numeric_limits<RealType>::quiet_NaN(), 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(0, 0))));
 706          BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(1, 0))));
 707          BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(2, 0))));
 708          BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(static_cast<RealType>(2.0001L), 0))));
 709          BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(3, 0))));
 710          BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(4, 0))));
 711
 712          // Kurtosis excess
 713          BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(std::numeric_limits<RealType>::quiet_NaN(), 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(0, 0))));
 716          BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(1, 0))));
 717          BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(2, 0))));
 718          BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(static_cast<RealType>(2.0001L), 0))));
 719          BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(3, 0))));
 720          BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(4, 0))));
 721       } // has_quiet_NaN
 722       BOOST_CHECK(boost::math::isfinite(mean(ignore_error_non_central_t(1 + std::numeric_limits<RealType>::epsilon(), 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(variance(ignore_error_non_central_t(static_cast<RealType>(2.0001L), 0))));
 725       BOOST_CHECK(boost::math::isfinite(variance(ignore_error_non_central_t(2 + 2 * std::numeric_limits<RealType>::epsilon(), 0))));
 726       BOOST_CHECK(boost::math::isfinite(skewness(ignore_error_non_central_t(3 + 3 * std::numeric_limits<RealType>::epsilon(), 0))));
 727       BOOST_CHECK(boost::math::isfinite(kurtosis(ignore_error_non_central_t(4 + 4 * std::numeric_limits<RealType>::epsilon(), 0))));
 728       BOOST_CHECK(boost::math::isfinite(kurtosis(ignore_error_non_central_t(static_cast<RealType>(4.0001L), 0))));
 729
 730       // check_out_of_range<non_central_t_distribution<RealType> >(1, 0); // Fails one check because allows df = infinity.
 731       check_support<non_central_t_distribution<RealType> >(non_central_t_distribution<RealType>(1, 0));
 732    } // ordinary floats.
 733 } // template <class RealType> void test_ignore_policy(RealType)
 734