math/test/test_nc_t.cpp
2019-08-10 08:50:12 -04:00

503 lines
19 KiB
C++

// test_nc_t.cpp
// Copyright John Maddock 2008, 2012.
// Copyright Paul A. Bristow 2012.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <pch.hpp> // Need to include lib/math/test in path.
#ifdef _MSC_VER
#pragma warning (disable:4127 4512)
#endif
#if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
# define TEST_FLOAT
# define TEST_DOUBLE
# define TEST_LDOUBLE
# define TEST_REAL_CONCEPT
#endif
#include <boost/math/tools/test.hpp>
#include <boost/math/concepts/real_concept.hpp> // for real_concept
#include <boost/math/distributions/non_central_t.hpp> // for chi_squared_distribution.
#include <boost/math/distributions/normal.hpp> // for normal distribution (for comparison).
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp> // for test_main
#include <boost/test/results_collector.hpp>
#include <boost/test/unit_test.hpp>
#include <boost/test/tools/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE
#include "functor.hpp"
#include "handle_test_result.hpp"
#include "table_type.hpp"
#include "test_nc_t.hpp"
#include <iostream>
#include <iomanip>
using std::cout;
using std::endl;
#include <limits>
using std::numeric_limits;
void expected_results()
{
//
// Define the max and mean errors expected for
// various compilers and platforms.
//
const char* largest_type;
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
if(boost::math::policies::digits<double, boost::math::policies::policy<> >() == boost::math::policies::digits<long double, boost::math::policies::policy<> >())
{
largest_type = "(long\\s+)?double|real_concept";
}
else
{
largest_type = "long double|real_concept";
}
#else
largest_type = "(long\\s+)?double|real_concept";
#endif
//
// Catch all cases come last:
//
if(std::numeric_limits<long double>::digits > 54)
{
add_expected_result(
"[^|]*", // compiler
"[^|]*", // stdlib
"[^|]*", // platform
largest_type, // test type(s)
"[^|]*large[^|]*", // test data group
"[^|]*", 2000000, 200000); // test function
add_expected_result(
"[^|]*", // compiler
"[^|]*", // stdlib
"[^|]*", // platform
"double", // test type(s)
"[^|]*large[^|]*", // test data group
"[^|]*", 500, 100); // test function
}
add_expected_result(
"[^|]*", // compiler
"[^|]*", // stdlib
"[^|]*", // platform
"real_concept", // test type(s)
"[^|]*", // test data group
"[^|]*", 300000, 100000); // test function
add_expected_result(
"[^|]*", // compiler
"[^|]*", // stdlib
"[^|]*", // platform
largest_type, // test type(s)
"[^|]*large[^|]*", // test data group
"[^|]*", 1500, 300); // test function
add_expected_result(
"[^|]*", // compiler
"[^|]*", // stdlib
"[^|]*", // platform
largest_type, // test type(s)
"[^|]*small[^|]*", // test data group
"[^|]*", 400, 100); // test function
add_expected_result(
"[^|]*", // compiler
"[^|]*", // stdlib
".*Solaris.*", // platform
largest_type, // test type(s)
"[^|]*", // test data group
"[^|]*", 400, 100); // test function
add_expected_result(
"[^|]*", // compiler
"[^|]*", // stdlib
"[^|]*", // platform
largest_type, // test type(s)
"[^|]*", // test data group
"[^|]*", 250, 50); // test function
//
// Finish off by printing out the compiler/stdlib/platform names,
// we do this to make it easier to mark up expected error rates.
//
std::cout << "Tests run with " << BOOST_COMPILER << ", "
<< BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl;
}
BOOST_AUTO_TEST_CASE( test_main )
{
BOOST_MATH_CONTROL_FP;
// Basic sanity-check spot values.
expected_results();
// (Parameter value, arbitrarily zero, only communicates the floating point type).
#ifdef TEST_FLOAT
test_spots(0.0F); // Test float.
#endif
#ifdef TEST_DOUBLE
test_spots(0.0); // Test double.
#endif
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
#ifdef TEST_LDOUBLE
test_spots(0.0L); // Test long double.
#endif
#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
#ifdef TEST_REAL_CONCEPT
test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
#endif
#endif
#endif
#ifdef TEST_FLOAT
test_accuracy(0.0F, "float"); // Test float.
test_big_df(0.F); // float
#endif
#ifdef TEST_DOUBLE
test_accuracy(0.0, "double"); // Test double.
test_big_df(0.); // double
test_ignore_policy(0.0);
#endif
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
#ifdef TEST_LDOUBLE
test_accuracy(0.0L, "long double"); // Test long double.
#endif
#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
#ifdef TEST_REAL_CONCEPT
test_accuracy(boost::math::concepts::real_concept(0.), "real_concept"); // Test real concept.
#endif
#endif
#endif
/* */
} // BOOST_AUTO_TEST_CASE( test_main )
/*
Output:
Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Debug\test_nc_t.exe"
Running 1 test case...
Tests run with Microsoft Visual C++ version 10.0, Dinkumware standard library version 520, Win32
Tolerance = 0.000596046%.
Tolerance = 5e-010%.
Tolerance = 5e-010%.
Tolerance = 1e-008%.
Testing: Non Central T
CDF<float> Max = 0 RMS Mean=0
CCDF<float> Max = 0 RMS Mean=0
Testing: float quantile sanity check, with tests Non Central T
Testing: Non Central T (small non-centrality)
CDF<float> Max = 0 RMS Mean=0
CCDF<float> Max = 0 RMS Mean=0
Testing: float quantile sanity check, with tests Non Central T (small non-centrality)
Testing: Non Central T (large parameters)
CDF<float> Max = 0 RMS Mean=0
CCDF<float> Max = 0 RMS Mean=0
Testing: float quantile sanity check, with tests Non Central T (large parameters)
Testing: Non Central T
CDF<double> Max = 137.7 RMS Mean=31.5
worst case at row: 181
{ 188.01481628417969, -282.022216796875, -298.02532958984375, 0.1552789395983287, 0.84472106040167128 }
CCDF<double> Max = 150.4 RMS Mean=32.32
worst case at row: 184
{ 191.43339538574219, 765.73358154296875, 820.14422607421875, 0.89943076553533785, 0.10056923446466212 }
Testing: double quantile sanity check, with tests Non Central T
Testing: Non Central T (small non-centrality)
CDF<double> Max = 3.605 RMS Mean=1.031
worst case at row: 42
{ 7376104448, 7.3761043495323975e-007, -1.3614851236343384, 0.086680099352107118, 0.91331990064789292 }
CCDF<double> Max = 5.207 RMS Mean=1.432
worst case at row: 38
{ 1524088576, 1.5240885886669275e-007, 1.3784774541854858, 0.91597201432644526, 0.084027985673554725 }
Testing: double quantile sanity check, with tests Non Central T (small non-centrality)
Testing: Non Central T (large parameters)
CDF<double> Max = 286.4 RMS Mean=62.79
worst case at row: 24
{ 1.3091821180254421e+019, 1309.18212890625, 1308.01171875, 0.12091797523015677, 0.87908202476984321 }
CCDF<double> Max = 226.9 RMS Mean=50.41
worst case at row: 23
{ 7.9217674231144776e+018, 792.1767578125, 793.54827880859375, 0.91489369852628, 0.085106301473719961 }
Testing: double quantile sanity check, with tests Non Central T (large parameters)
Testing: Non Central T
CDF<long double> Max = 137.7 RMS Mean=31.5
worst case at row: 181
{ 188.01481628417969, -282.022216796875, -298.02532958984375, 0.1552789395983287, 0.84472106040167128 }
CCDF<long double> Max = 150.4 RMS Mean=32.32
worst case at row: 184
{ 191.43339538574219, 765.73358154296875, 820.14422607421875, 0.89943076553533785, 0.10056923446466212 }
Testing: long double quantile sanity check, with tests Non Central T
Testing: Non Central T (small non-centrality)
CDF<long double> Max = 3.605 RMS Mean=1.031
worst case at row: 42
{ 7376104448, 7.3761043495323975e-007, -1.3614851236343384, 0.086680099352107118, 0.91331990064789292 }
CCDF<long double> Max = 5.207 RMS Mean=1.432
worst case at row: 38
{ 1524088576, 1.5240885886669275e-007, 1.3784774541854858, 0.91597201432644526, 0.084027985673554725 }
Testing: long double quantile sanity check, with tests Non Central T (small non-centrality)
Testing: Non Central T (large parameters)
CDF<long double> Max = 286.4 RMS Mean=62.79
worst case at row: 24
{ 1.3091821180254421e+019, 1309.18212890625, 1308.01171875, 0.12091797523015677, 0.87908202476984321 }
CCDF<long double> Max = 226.9 RMS Mean=50.41
worst case at row: 23
{ 7.9217674231144776e+018, 792.1767578125, 793.54827880859375, 0.91489369852628, 0.085106301473719961 }
Testing: long double quantile sanity check, with tests Non Central T (large parameters)
Testing: Non Central T
CDF<real_concept> Max = 2.816e+005 RMS Mean=2.029e+004
worst case at row: 185
{ 191.50137329101562, -957.5068359375, -1035.4078369140625, 0.072545502958829097, 0.92745449704117089 }
CCDF<real_concept> Max = 1.304e+005 RMS Mean=1.529e+004
worst case at row: 184
{ 191.43339538574219, 765.73358154296875, 820.14422607421875, 0.89943076553533785, 0.10056923446466212 }
cdf(n10, 11) = 0.84134471416473389 0.15865525603294373
cdf(n10, 9) = 0.15865525603294373 0.84134471416473389
cdf(maxdf10, 11) = 0.84134477376937866 0.15865525603294373
cdf(infdf10, 11) = 0.84134477376937866 0.15865525603294373
cdf(n10, 11) = 0.84134474606854293 0.15865525393145707
cdf(n10, 9) = 0.15865525393145707 0.84134474606854293
cdf(maxdf10, 11) = 0.84134474606854293 0.15865525393145707
cdf(infdf10, 11) = 0.84134474606854293 0.15865525393145707
*** No errors detected
Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Debug\test_nc_t.exe"
Running 1 test case...
Tests run with Microsoft Visual C++ version 10.0, Dinkumware standard library version 520, Win32
Tolerance = 0.000596046%.
Tolerance = 5e-010%.
Tolerance = 5e-010%.
Tolerance = 1e-008%.
Testing: Non Central T
CDF<float> Max = 0 RMS Mean=0
CCDF<float> Max = 0 RMS Mean=0
Testing: float quantile sanity check, with tests Non Central T
Testing: Non Central T (small non-centrality)
CDF<float> Max = 0 RMS Mean=0
CCDF<float> Max = 0 RMS Mean=0
Testing: float quantile sanity check, with tests Non Central T (small non-centrality)
Testing: Non Central T (large parameters)
CDF<float> Max = 0 RMS Mean=0
CCDF<float> Max = 0 RMS Mean=0
Testing: float quantile sanity check, with tests Non Central T (large parameters)
Testing: Non Central T
CDF<double> Max = 137.7 RMS Mean=31.5
worst case at row: 181
{ 188.01481628417969, -282.022216796875, -298.02532958984375, 0.1552789395983287, 0.84472106040167128 }
CCDF<double> Max = 150.4 RMS Mean=32.32
worst case at row: 184
{ 191.43339538574219, 765.73358154296875, 820.14422607421875, 0.89943076553533785, 0.10056923446466212 }
Testing: double quantile sanity check, with tests Non Central T
Testing: Non Central T (small non-centrality)
CDF<double> Max = 3.605 RMS Mean=1.031
worst case at row: 42
{ 7376104448, 7.3761043495323975e-007, -1.3614851236343384, 0.086680099352107118, 0.91331990064789292 }
CCDF<double> Max = 5.207 RMS Mean=1.432
worst case at row: 38
{ 1524088576, 1.5240885886669275e-007, 1.3784774541854858, 0.91597201432644526, 0.084027985673554725 }
Testing: double quantile sanity check, with tests Non Central T (small non-centrality)
Testing: Non Central T (large parameters)
CDF<double> Max = 286.4 RMS Mean=62.79
worst case at row: 24
{ 1.3091821180254421e+019, 1309.18212890625, 1308.01171875, 0.12091797523015677, 0.87908202476984321 }
CCDF<double> Max = 226.9 RMS Mean=50.41
worst case at row: 23
{ 7.9217674231144776e+018, 792.1767578125, 793.54827880859375, 0.91489369852628, 0.085106301473719961 }
Testing: double quantile sanity check, with tests Non Central T (large parameters)
Testing: Non Central T
CDF<long double> Max = 137.7 RMS Mean=31.5
worst case at row: 181
{ 188.01481628417969, -282.022216796875, -298.02532958984375, 0.1552789395983287, 0.84472106040167128 }
CCDF<long double> Max = 150.4 RMS Mean=32.32
worst case at row: 184
{ 191.43339538574219, 765.73358154296875, 820.14422607421875, 0.89943076553533785, 0.10056923446466212 }
Testing: long double quantile sanity check, with tests Non Central T
Testing: Non Central T (small non-centrality)
CDF<long double> Max = 3.605 RMS Mean=1.031
worst case at row: 42
{ 7376104448, 7.3761043495323975e-007, -1.3614851236343384, 0.086680099352107118, 0.91331990064789292 }
CCDF<long double> Max = 5.207 RMS Mean=1.432
worst case at row: 38
{ 1524088576, 1.5240885886669275e-007, 1.3784774541854858, 0.91597201432644526, 0.084027985673554725 }
Testing: long double quantile sanity check, with tests Non Central T (small non-centrality)
Testing: Non Central T (large parameters)
CDF<long double> Max = 286.4 RMS Mean=62.79
worst case at row: 24
{ 1.3091821180254421e+019, 1309.18212890625, 1308.01171875, 0.12091797523015677, 0.87908202476984321 }
CCDF<long double> Max = 226.9 RMS Mean=50.41
worst case at row: 23
{ 7.9217674231144776e+018, 792.1767578125, 793.54827880859375, 0.91489369852628, 0.085106301473719961 }
Testing: long double quantile sanity check, with tests Non Central T (large parameters)
Testing: Non Central T
CDF<real_concept> Max = 2.816e+005 RMS Mean=2.029e+004
worst case at row: 185
{ 191.50137329101562, -957.5068359375, -1035.4078369140625, 0.072545502958829097, 0.92745449704117089 }
CCDF<real_concept> Max = 1.304e+005 RMS Mean=1.529e+004
worst case at row: 184
{ 191.43339538574219, 765.73358154296875, 820.14422607421875, 0.89943076553533785, 0.10056923446466212 }
*** No errors detected
*/
/*
Temporary stuff from student's t version.
// Calculate 1 / eps, the point where student's t should change to normal distribution.
RealType limit = 1 / boost::math::tools::epsilon<RealType>();
using namespace boost::math::policies;
typedef policy<digits10<17> > accurate_policy; // 17 = max_digits10 where available.
limit = 1 / policies::get_epsilon<RealType, accurate_policy>();
BOOST_CHECK_CLOSE_FRACTION(limit, static_cast<RealType>(1) / std::numeric_limits<RealType>::epsilon(), tolerance);
// Default policy to get full accuracy.
// std::cout << "Switch over to normal if df > " << limit << std::endl;
// float Switch over to normal if df > 8.38861e+006
// double Switch over to normal if df > 4.5036e+015
// Can't test real_concept - doesn't converge.
boost::math::normal_distribution<RealType> n01(0, 1); //
boost::math::normal_distribution<RealType> n10(10, 1); //
non_central_t_distribution<RealType> nct(boost::math::tools::max_value<RealType>(), 0); // Well over the switchover point,
non_central_t_distribution<RealType> nct2(limit /5, 0); // Just below the switchover point,
non_central_t_distribution<RealType> nct3(limit /100, 0); // Well below the switchover point,
non_central_t_distribution<RealType> nct4(limit, 10); // Well below the switchover point, and 10 non-centrality.
// PDF
BOOST_CHECK_CLOSE_FRACTION(pdf(nct, 0), pdf(n01, 0.), tolerance); // normal and non-central t should be nearly equal.
BOOST_CHECK_CLOSE_FRACTION(pdf(nct2, 0), pdf(n01, 0.), tolerance); // should be very close to normal.
BOOST_CHECK_CLOSE_FRACTION(pdf(nct3, 0), pdf(n01, 0.), tolerance * 10); // should be close to normal.
// BOOST_CHECK_CLOSE_FRACTION(pdf(nct4, 10), pdf(n10, 0.), tolerance * 100); // should be fairly close to normal tolerance.
RealType delta = 10; // non-centrality.
RealType nu = static_cast<RealType>(limit); // df
boost::math::normal_distribution<RealType> nl(delta, 1); // Normal distribution that nct tends to for big df.
non_central_t_distribution<RealType> nct5(nu, delta); //
RealType x = delta;
// BOOST_CHECK_CLOSE_FRACTION(pdf(nct5, x), pdf(nl, x), tolerance * 10 ); // nu = 1e15
// BOOST_CHECK_CLOSE_FRACTION(pdf(nct5, x), pdf(nl, x), tolerance * 1000 ); // nu = 1e14
// BOOST_CHECK_CLOSE_FRACTION(pdf(nct5, x), pdf(nl, x), tolerance * 10000 ); // nu = 1e13
// BOOST_CHECK_CLOSE_FRACTION(pdf(nct5, x), pdf(nl, x), tolerance * 100000 ); // nu = 1e12
BOOST_CHECK_CLOSE_FRACTION(pdf(nct5, x), pdf(nl, x), tolerance * 5 ); // nu = 1/eps
// Increasing the non-centrality delta increases the difference too because increases asymmetry.
// For example, with non-centrality = 100, need tolerance * 500
// CDF
BOOST_CHECK_CLOSE_FRACTION(cdf(nct, 0), cdf(n01, 0.), tolerance); // should be exactly equal.
BOOST_CHECK_CLOSE_FRACTION(cdf(nct2, 0), cdf(n01, 0.), tolerance); // should be very close to normal.
BOOST_CHECK_CLOSE_FRACTION(cdf(complement(n10, 11)), 1 - cdf(n10, 11), tolerance); //
// cdf(n10, 10) = 0.841345 0.158655
BOOST_CHECK_CLOSE_FRACTION(cdf(complement(n10, 9)), 1 - cdf(n10, 9), tolerance); //
std::cout.precision(17);
std::cout << "cdf(n10, 11) = " << cdf(n10, 11) << ' ' << cdf(complement(n10, 11)) << endl;
std::cout << "cdf(n10, 9) = " << cdf(n10, 9) << ' ' << cdf(complement(n10, 9)) << endl;
std::cout << std::numeric_limits<double>::max_digits10 << std::endl;
std::cout.precision(17);
using boost::math::tools::max_value;
double eps = std::numeric_limits<double>::epsilon();
// Use policies so that if policy requests lower precision,
// then get the normal distribution approximation earlier.
//limit = static_cast<double>(1) / limit; // 1/eps
double delta = 1e2;
double df =
delta / (4 * eps);
std::cout << df << std::endl; // df = 1.125899906842624e+018
{
boost::math::non_central_t_distribution<double> dist(df, delta);
std::cout <<"mean " << mean(dist) << std::endl; // mean 1000
std::cout <<"variance " << variance(dist) << std::endl; // variance 1
std::cout <<"skewness " << skewness(dist) << std::endl; // skewness 8.8817841970012523e-010
std::cout <<"kurtosis_excess " << kurtosis_excess(dist) << std::endl; // kurtosis_excess 3.0001220703125
//1.125899906842624e+017
//mean 100
//variance 1
//skewness 8.8817841970012523e-012
//kurtosis_excess 3
}
*/