libs/math/test/test_nc_t.cpp

   1 // test_nc_t.cpp
   2
   3 // Copyright John Maddock 2008, 2012.
   4 // Copyright Paul A. Bristow 2012.
   5
   6 // Use, modification and distribution are subject to the
   7 // Boost Software License, Version 1.0.
   8 // (See accompanying file LICENSE_1_0.txt
   9 // or copy at http://www.boost.org/LICENSE_1_0.txt)
  10
  11 #include <pch.hpp> // Need to include lib/math/test in path.
  12
  13 #ifdef _MSC_VER
  14 #pragma warning (disable:4127 4512)
  15 #endif
  16
  17 #if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
  18 #  define TEST_FLOAT
  19 #  define TEST_DOUBLE
  20 #  define TEST_LDOUBLE
  21 #  define TEST_REAL_CONCEPT
  22 #endif
  23
  24 #include <boost/math/tools/test.hpp>
  25 #include <boost/math/concepts/real_concept.hpp> // for real_concept
  26 #include <boost/math/distributions/non_central_t.hpp> // for chi_squared_distribution.
  27 #include <boost/math/distributions/normal.hpp> // for normal distribution (for comparison).
  28
  29 #define BOOST_TEST_MAIN
  30 #include <boost/test/unit_test.hpp> // for test_main
  31 #include <boost/test/results_collector.hpp>
  32 #include <boost/test/unit_test.hpp>
  33 #include <boost/test/tools/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE
  34
  35 #include "functor.hpp"
  36 #include "handle_test_result.hpp"
  37 #include "table_type.hpp"
  38 #include "test_nc_t.hpp"
  39
  40 #include <iostream>
  41 #include <iomanip>
  42 using std::cout;
  43 using std::endl;
  44 #include <limits>
  45 using std::numeric_limits;
  46
  47
  48 void expected_results()
  49 {
  50    //
  51    // Define the max and mean errors expected for
  52    // various compilers and platforms.
  53    //
  54    const char* largest_type;
  55 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  56    if(boost::math::policies::digits<double, boost::math::policies::policy<> >() == boost::math::policies::digits<long double, boost::math::policies::policy<> >())
  57    {
  58       largest_type = "(long\\s+)?double|real_concept";
  59    }
  60    else
  61    {
  62       largest_type = "long double|real_concept";
  63    }
  64 #else
  65    largest_type = "(long\\s+)?double|real_concept";
  66 #endif
  67
  68    //
  69    // Catch all cases come last:
  70    //
  71    if(std::numeric_limits<long double>::digits > 54)
  72    {
  73       add_expected_result(
  74          "[^|]*",                          // compiler
  75          "[^|]*",                          // stdlib
  76          "[^|]*",                          // platform
  77          largest_type,                     // test type(s)
  78          "[^|]*large[^|]*",                // test data group
  79          "[^|]*", 2000000, 200000);        // test function
  80       add_expected_result(
  81          "[^|]*",                          // compiler
  82          "[^|]*",                          // stdlib
  83          "[^|]*",                          // platform
  84          "double",                         // test type(s)
  85          "[^|]*large[^|]*",                // test data group
  86          "[^|]*", 500, 100);               // test function
  87    }
  88    add_expected_result(
  89       "[^|]*",                          // compiler
  90       "[^|]*",                          // stdlib
  91       "[^|]*",                          // platform
  92       "real_concept",                   // test type(s)
  93       "[^|]*",                          // test data group
  94       "[^|]*", 300000, 100000);                // test function
  95    add_expected_result(
  96       "[^|]*",                          // compiler
  97       "[^|]*",                          // stdlib
  98       "[^|]*",                          // platform
  99       largest_type,                     // test type(s)
 100       "[^|]*large[^|]*",                // test data group
 101       "[^|]*", 1500, 300);              // test function
 102    add_expected_result(
 103       "[^|]*",                          // compiler
 104       "[^|]*",                          // stdlib
 105       "[^|]*",                          // platform
 106       largest_type,                     // test type(s)
 107       "[^|]*small[^|]*",                // test data group
 108       "[^|]*", 400, 100);              // test function
 109    add_expected_result(
 110       "[^|]*",                          // compiler
 111       "[^|]*",                          // stdlib
 112       ".*Solaris.*",                    // platform
 113       largest_type,                     // test type(s)
 114       "[^|]*",                          // test data group
 115       "[^|]*", 400, 100);               // test function
 116    add_expected_result(
 117       "[^|]*",                          // compiler
 118       "[^|]*",                          // stdlib
 119       "[^|]*",                          // platform
 120       largest_type,                     // test type(s)
 121       "[^|]*",                          // test data group
 122       "[^|]*", 250, 50);                // test function
 123
 124    //
 125    // Finish off by printing out the compiler/stdlib/platform names,
 126    // we do this to make it easier to mark up expected error rates.
 127    //
 128    std::cout << "Tests run with " << BOOST_COMPILER << ", "
 129       << BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl;
 130 }
 131
 132
 133 BOOST_AUTO_TEST_CASE( test_main )
 134 {
 135   BOOST_MATH_CONTROL_FP;
 136    // Basic sanity-check spot values.
 137    expected_results();
 138
 139    // (Parameter value, arbitrarily zero, only communicates the floating point type).
 140 #ifdef TEST_FLOAT
 141    test_spots(0.0F); // Test float.
 142 #endif
 143 #ifdef TEST_DOUBLE
 144    test_spots(0.0); // Test double.
 145 #endif
 146 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 147 #ifdef TEST_LDOUBLE
 148    test_spots(0.0L); // Test long double.
 149 #endif
 150 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
 151 #ifdef TEST_REAL_CONCEPT
 152    test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
 153 #endif
 154 #endif
 155 #endif
 156
 157 #ifdef TEST_FLOAT
 158    test_accuracy(0.0F, "float"); // Test float.
 159    test_big_df(0.F); // float
 160 #endif
 161 #ifdef TEST_DOUBLE
 162    test_accuracy(0.0, "double"); // Test double.
 163    test_big_df(0.); // double
 164    test_ignore_policy(0.0);
 165 #endif
 166 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 167 #ifdef TEST_LDOUBLE
 168    test_accuracy(0.0L, "long double"); // Test long double.
 169 #endif
 170 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
 171 #ifdef TEST_REAL_CONCEPT
 172    test_accuracy(boost::math::concepts::real_concept(0.), "real_concept"); // Test real concept.
 173 #endif
 174 #endif
 175 #endif
 176   /* */
 177
 178
 179 } // BOOST_AUTO_TEST_CASE( test_main )
 180
 181 /*
 182
 183 Output:
 184
 185   Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Debug\test_nc_t.exe"
 186   Running 1 test case...
 187   Tests run with Microsoft Visual C++ version 10.0, Dinkumware standard library version 520, Win32
 188   Tolerance = 0.000596046%.
 189   Tolerance = 5e-010%.
 190   Tolerance = 5e-010%.
 191   Tolerance = 1e-008%.
 192   Testing: Non Central T
 193   CDF<float> Max = 0 RMS Mean=0
 194
 195   CCDF<float> Max = 0 RMS Mean=0
 196
 197
 198   Testing: float quantile sanity check, with tests Non Central T
 199   Testing: Non Central T (small non-centrality)
 200   CDF<float> Max = 0 RMS Mean=0
 201
 202   CCDF<float> Max = 0 RMS Mean=0
 203
 204
 205   Testing: float quantile sanity check, with tests Non Central T (small non-centrality)
 206   Testing: Non Central T (large parameters)
 207   CDF<float> Max = 0 RMS Mean=0
 208
 209   CCDF<float> Max = 0 RMS Mean=0
 210
 211
 212   Testing: float quantile sanity check, with tests Non Central T (large parameters)
 213   Testing: Non Central T
 214   CDF<double> Max = 137.7 RMS Mean=31.5
 215       worst case at row: 181
 216       { 188.01481628417969, -282.022216796875, -298.02532958984375, 0.1552789395983287, 0.84472106040167128 }
 217
 218   CCDF<double> Max = 150.4 RMS Mean=32.32
 219       worst case at row: 184
 220       { 191.43339538574219, 765.73358154296875, 820.14422607421875, 0.89943076553533785, 0.10056923446466212 }
 221
 222
 223   Testing: double quantile sanity check, with tests Non Central T
 224   Testing: Non Central T (small non-centrality)
 225   CDF<double> Max = 3.605 RMS Mean=1.031
 226       worst case at row: 42
 227       { 7376104448, 7.3761043495323975e-007, -1.3614851236343384, 0.086680099352107118, 0.91331990064789292 }
 228
 229   CCDF<double> Max = 5.207 RMS Mean=1.432
 230       worst case at row: 38
 231       { 1524088576, 1.5240885886669275e-007, 1.3784774541854858, 0.91597201432644526, 0.084027985673554725 }
 232
 233
 234   Testing: double quantile sanity check, with tests Non Central T (small non-centrality)
 235   Testing: Non Central T (large parameters)
 236   CDF<double> Max = 286.4 RMS Mean=62.79
 237       worst case at row: 24
 238       { 1.3091821180254421e+019, 1309.18212890625, 1308.01171875, 0.12091797523015677, 0.87908202476984321 }
 239
 240   CCDF<double> Max = 226.9 RMS Mean=50.41
 241       worst case at row: 23
 242       { 7.9217674231144776e+018, 792.1767578125, 793.54827880859375, 0.91489369852628, 0.085106301473719961 }
 243
 244
 245   Testing: double quantile sanity check, with tests Non Central T (large parameters)
 246   Testing: Non Central T
 247   CDF<long double> Max = 137.7 RMS Mean=31.5
 248       worst case at row: 181
 249       { 188.01481628417969, -282.022216796875, -298.02532958984375, 0.1552789395983287, 0.84472106040167128 }
 250
 251   CCDF<long double> Max = 150.4 RMS Mean=32.32
 252       worst case at row: 184
 253       { 191.43339538574219, 765.73358154296875, 820.14422607421875, 0.89943076553533785, 0.10056923446466212 }
 254
 255
 256   Testing: long double quantile sanity check, with tests Non Central T
 257   Testing: Non Central T (small non-centrality)
 258   CDF<long double> Max = 3.605 RMS Mean=1.031
 259       worst case at row: 42
 260       { 7376104448, 7.3761043495323975e-007, -1.3614851236343384, 0.086680099352107118, 0.91331990064789292 }
 261
 262   CCDF<long double> Max = 5.207 RMS Mean=1.432
 263       worst case at row: 38
 264       { 1524088576, 1.5240885886669275e-007, 1.3784774541854858, 0.91597201432644526, 0.084027985673554725 }
 265
 266
 267   Testing: long double quantile sanity check, with tests Non Central T (small non-centrality)
 268   Testing: Non Central T (large parameters)
 269   CDF<long double> Max = 286.4 RMS Mean=62.79
 270       worst case at row: 24
 271       { 1.3091821180254421e+019, 1309.18212890625, 1308.01171875, 0.12091797523015677, 0.87908202476984321 }
 272
 273   CCDF<long double> Max = 226.9 RMS Mean=50.41
 274       worst case at row: 23
 275       { 7.9217674231144776e+018, 792.1767578125, 793.54827880859375, 0.91489369852628, 0.085106301473719961 }
 276
 277
 278   Testing: long double quantile sanity check, with tests Non Central T (large parameters)
 279   Testing: Non Central T
 280   CDF<real_concept> Max = 2.816e+005 RMS Mean=2.029e+004
 281       worst case at row: 185
 282       { 191.50137329101562, -957.5068359375, -1035.4078369140625, 0.072545502958829097, 0.92745449704117089 }
 283
 284   CCDF<real_concept> Max = 1.304e+005 RMS Mean=1.529e+004
 285       worst case at row: 184
 286       { 191.43339538574219, 765.73358154296875, 820.14422607421875, 0.89943076553533785, 0.10056923446466212 }
 287
 288
 289   cdf(n10, 11)  = 0.84134471416473389 0.15865525603294373
 290   cdf(n10, 9)  = 0.15865525603294373 0.84134471416473389
 291   cdf(maxdf10, 11)  = 0.84134477376937866 0.15865525603294373
 292   cdf(infdf10, 11)  = 0.84134477376937866 0.15865525603294373
 293   cdf(n10, 11)  = 0.84134474606854293 0.15865525393145707
 294   cdf(n10, 9)  = 0.15865525393145707 0.84134474606854293
 295   cdf(maxdf10, 11)  = 0.84134474606854293 0.15865525393145707
 296   cdf(infdf10, 11)  = 0.84134474606854293 0.15865525393145707
 297
 298   *** No errors detected
 299
 300     Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Debug\test_nc_t.exe"
 301   Running 1 test case...
 302   Tests run with Microsoft Visual C++ version 10.0, Dinkumware standard library version 520, Win32
 303   Tolerance = 0.000596046%.
 304   Tolerance = 5e-010%.
 305   Tolerance = 5e-010%.
 306   Tolerance = 1e-008%.
 307   Testing: Non Central T
 308   CDF<float> Max = 0 RMS Mean=0
 309
 310   CCDF<float> Max = 0 RMS Mean=0
 311
 312
 313   Testing: float quantile sanity check, with tests Non Central T
 314   Testing: Non Central T (small non-centrality)
 315   CDF<float> Max = 0 RMS Mean=0
 316
 317   CCDF<float> Max = 0 RMS Mean=0
 318
 319
 320   Testing: float quantile sanity check, with tests Non Central T (small non-centrality)
 321   Testing: Non Central T (large parameters)
 322   CDF<float> Max = 0 RMS Mean=0
 323
 324   CCDF<float> Max = 0 RMS Mean=0
 325
 326
 327   Testing: float quantile sanity check, with tests Non Central T (large parameters)
 328   Testing: Non Central T
 329   CDF<double> Max = 137.7 RMS Mean=31.5
 330       worst case at row: 181
 331       { 188.01481628417969, -282.022216796875, -298.02532958984375, 0.1552789395983287, 0.84472106040167128 }
 332
 333   CCDF<double> Max = 150.4 RMS Mean=32.32
 334       worst case at row: 184
 335       { 191.43339538574219, 765.73358154296875, 820.14422607421875, 0.89943076553533785, 0.10056923446466212 }
 336
 337
 338   Testing: double quantile sanity check, with tests Non Central T
 339   Testing: Non Central T (small non-centrality)
 340   CDF<double> Max = 3.605 RMS Mean=1.031
 341       worst case at row: 42
 342       { 7376104448, 7.3761043495323975e-007, -1.3614851236343384, 0.086680099352107118, 0.91331990064789292 }
 343
 344   CCDF<double> Max = 5.207 RMS Mean=1.432
 345       worst case at row: 38
 346       { 1524088576, 1.5240885886669275e-007, 1.3784774541854858, 0.91597201432644526, 0.084027985673554725 }
 347
 348
 349   Testing: double quantile sanity check, with tests Non Central T (small non-centrality)
 350   Testing: Non Central T (large parameters)
 351   CDF<double> Max = 286.4 RMS Mean=62.79
 352       worst case at row: 24
 353       { 1.3091821180254421e+019, 1309.18212890625, 1308.01171875, 0.12091797523015677, 0.87908202476984321 }
 354
 355   CCDF<double> Max = 226.9 RMS Mean=50.41
 356       worst case at row: 23
 357       { 7.9217674231144776e+018, 792.1767578125, 793.54827880859375, 0.91489369852628, 0.085106301473719961 }
 358
 359
 360   Testing: double quantile sanity check, with tests Non Central T (large parameters)
 361   Testing: Non Central T
 362   CDF<long double> Max = 137.7 RMS Mean=31.5
 363       worst case at row: 181
 364       { 188.01481628417969, -282.022216796875, -298.02532958984375, 0.1552789395983287, 0.84472106040167128 }
 365
 366   CCDF<long double> Max = 150.4 RMS Mean=32.32
 367       worst case at row: 184
 368       { 191.43339538574219, 765.73358154296875, 820.14422607421875, 0.89943076553533785, 0.10056923446466212 }
 369
 370
 371   Testing: long double quantile sanity check, with tests Non Central T
 372   Testing: Non Central T (small non-centrality)
 373   CDF<long double> Max = 3.605 RMS Mean=1.031
 374       worst case at row: 42
 375       { 7376104448, 7.3761043495323975e-007, -1.3614851236343384, 0.086680099352107118, 0.91331990064789292 }
 376
 377   CCDF<long double> Max = 5.207 RMS Mean=1.432
 378       worst case at row: 38
 379       { 1524088576, 1.5240885886669275e-007, 1.3784774541854858, 0.91597201432644526, 0.084027985673554725 }
 380
 381
 382   Testing: long double quantile sanity check, with tests Non Central T (small non-centrality)
 383   Testing: Non Central T (large parameters)
 384   CDF<long double> Max = 286.4 RMS Mean=62.79
 385       worst case at row: 24
 386       { 1.3091821180254421e+019, 1309.18212890625, 1308.01171875, 0.12091797523015677, 0.87908202476984321 }
 387
 388   CCDF<long double> Max = 226.9 RMS Mean=50.41
 389       worst case at row: 23
 390       { 7.9217674231144776e+018, 792.1767578125, 793.54827880859375, 0.91489369852628, 0.085106301473719961 }
 391
 392
 393   Testing: long double quantile sanity check, with tests Non Central T (large parameters)
 394   Testing: Non Central T
 395   CDF<real_concept> Max = 2.816e+005 RMS Mean=2.029e+004
 396       worst case at row: 185
 397       { 191.50137329101562, -957.5068359375, -1035.4078369140625, 0.072545502958829097, 0.92745449704117089 }
 398
 399   CCDF<real_concept> Max = 1.304e+005 RMS Mean=1.529e+004
 400       worst case at row: 184
 401       { 191.43339538574219, 765.73358154296875, 820.14422607421875, 0.89943076553533785, 0.10056923446466212 }
 402
 403
 404
 405   *** No errors detected
 406
 407
 408 */
 409
 410
 411
 412 /*
 413
 414 Temporary stuff from student's t version.
 415
 416
 417    // Calculate 1 / eps, the point where student's t should change to normal distribution.
 418     RealType limit = 1 / boost::math::tools::epsilon<RealType>();
 419
 420     using namespace boost::math::policies;
 421     typedef policy<digits10<17> > accurate_policy; // 17 = max_digits10 where available.
 422     limit = 1 / policies::get_epsilon<RealType, accurate_policy>();
 423
 424     BOOST_CHECK_CLOSE_FRACTION(limit, static_cast<RealType>(1) / std::numeric_limits<RealType>::epsilon(), tolerance);
 425     // Default policy to get full accuracy.
 426     // std::cout << "Switch over to normal if df > " << limit << std::endl;
 427     // float Switch over to normal if df > 8.38861e+006
 428     // double Switch over to normal if df > 4.5036e+015
 429     // Can't test real_concept - doesn't converge.
 430
 431     boost::math::normal_distribution<RealType> n01(0, 1); //
 432     boost::math::normal_distribution<RealType> n10(10, 1); //
 433     non_central_t_distribution<RealType> nct(boost::math::tools::max_value<RealType>(), 0); // Well over the switchover point,
 434     non_central_t_distribution<RealType> nct2(limit /5, 0); // Just below the switchover point,
 435     non_central_t_distribution<RealType> nct3(limit /100, 0); // Well below the switchover point,
 436     non_central_t_distribution<RealType> nct4(limit, 10); // Well below the switchover point, and 10 non-centrality.
 437
 438     // PDF
 439     BOOST_CHECK_CLOSE_FRACTION(pdf(nct, 0), pdf(n01, 0.), tolerance); // normal and non-central t should be nearly equal.
 440     BOOST_CHECK_CLOSE_FRACTION(pdf(nct2, 0), pdf(n01, 0.), tolerance); // should be very close to normal.
 441     BOOST_CHECK_CLOSE_FRACTION(pdf(nct3, 0), pdf(n01, 0.), tolerance * 10); // should be close to normal.
 442  //   BOOST_CHECK_CLOSE_FRACTION(pdf(nct4, 10), pdf(n10, 0.), tolerance * 100); // should be fairly close to normal tolerance.
 443
 444     RealType delta = 10; // non-centrality.
 445     RealType nu = static_cast<RealType>(limit); // df
 446     boost::math::normal_distribution<RealType> nl(delta, 1); // Normal distribution that nct tends to for big df.
 447     non_central_t_distribution<RealType> nct5(nu, delta); //
 448     RealType x = delta;
 449   //  BOOST_CHECK_CLOSE_FRACTION(pdf(nct5, x), pdf(nl, x), tolerance * 10 ); // nu = 1e15
 450   //  BOOST_CHECK_CLOSE_FRACTION(pdf(nct5, x), pdf(nl, x), tolerance * 1000 ); // nu = 1e14
 451   //  BOOST_CHECK_CLOSE_FRACTION(pdf(nct5, x), pdf(nl, x), tolerance * 10000 ); // nu = 1e13
 452   //  BOOST_CHECK_CLOSE_FRACTION(pdf(nct5, x), pdf(nl, x), tolerance * 100000 ); // nu = 1e12
 453     BOOST_CHECK_CLOSE_FRACTION(pdf(nct5, x), pdf(nl, x), tolerance * 5  ); // nu = 1/eps
 454
 455   // Increasing the non-centrality delta increases the difference too because increases asymmetry.
 456   // For example, with non-centrality = 100, need tolerance * 500
 457
 458       // CDF
 459     BOOST_CHECK_CLOSE_FRACTION(cdf(nct, 0), cdf(n01, 0.), tolerance); // should be exactly equal.
 460     BOOST_CHECK_CLOSE_FRACTION(cdf(nct2, 0), cdf(n01, 0.), tolerance); // should be very close to normal.
 461
 462     BOOST_CHECK_CLOSE_FRACTION(cdf(complement(n10, 11)), 1 - cdf(n10, 11), tolerance); //
 463     //   cdf(n10, 10)  = 0.841345 0.158655
 464     BOOST_CHECK_CLOSE_FRACTION(cdf(complement(n10, 9)), 1 - cdf(n10, 9), tolerance); //
 465     std::cout.precision(17);
 466     std::cout  << "cdf(n10, 11)  = " << cdf(n10, 11) << ' ' << cdf(complement(n10, 11)) << endl;
 467     std::cout  << "cdf(n10, 9)  = " << cdf(n10, 9) << ' ' << cdf(complement(n10, 9)) << endl;
 468
 469   std::cout << std::numeric_limits<double>::max_digits10 << std::endl;
 470    std::cout.precision(17);
 471
 472    using boost::math::tools::max_value;
 473
 474    double eps = std::numeric_limits<double>::epsilon();
 475    // Use policies so that if policy requests lower precision,
 476    // then get the normal distribution approximation earlier.
 477    //limit = static_cast<double>(1) / limit; // 1/eps
 478    double delta = 1e2;
 479    double df =
 480    delta / (4 * eps);
 481
 482     std::cout << df << std::endl; // df = 1.125899906842624e+018
 483
 484    {
 485      boost::math::non_central_t_distribution<double> dist(df, delta);
 486
 487       std::cout <<"mean " << mean(dist) << std::endl; // mean 1000
 488       std::cout <<"variance " << variance(dist) << std::endl; // variance 1
 489       std::cout <<"skewness " << skewness(dist) << std::endl; //  skewness 8.8817841970012523e-010
 490       std::cout <<"kurtosis_excess " << kurtosis_excess(dist) << std::endl; // kurtosis_excess 3.0001220703125
 491   //1.125899906842624e+017
 492   //mean 100
 493   //variance 1
 494   //skewness 8.8817841970012523e-012
 495   //kurtosis_excess 3
 496
 497
 498    }
 499
 500
 501
 502   */