libs/math/test/test_owens_t.cpp

   1 // test_owens_t.cpp
   2
   3 // Copyright Paul A. Bristow 2012.
   4 // Copyright Benjamin Sobotta 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 // Tested using some 30 decimal digit accuracy values from:
  12 // Fast and accurate calculation of Owen's T-function
  13 // Mike Patefield, and David Tandy
  14 // Journal of Statistical Software, 5 (5), 1-25 (2000).
  15 // http://www.jstatsoft.org/v05/a05/paper  Table 3, page 15
  16 // Values of T(h,a) accurate to thirty figures were calculated using 128 bit arithmetic by
  17 // evaluating (9) with m = 48, the summation over k being continued until additional terms did
  18 // not alter the result. The resultant values Tacc(h,a) say, were validated by evaluating (8) with
  19 // m = 48 (i.e. 96 point Gaussian quadrature).
  20
  21 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
  22
  23 #ifdef _MSC_VER
  24 #  pragma warning (disable : 4127) // conditional expression is constant
  25 #  pragma warning (disable : 4305) // 'initializing' : truncation from 'double' to 'const float'
  26 // ?? TODO get rid of these warnings?
  27 #endif
  28
  29 #include <boost/math/concepts/real_concept.hpp> // for real_concept.
  30 using ::boost::math::concepts::real_concept;
  31
  32 #include <boost/math/special_functions/owens_t.hpp> // for owens_t function.
  33 using boost::math::owens_t;
  34 #include <boost/math/distributions/normal.hpp>
  35
  36 #define BOOST_TEST_MAIN
  37 #include <boost/test/unit_test.hpp>
  38 #include <boost/test/floating_point_comparison.hpp>
  39 #include <boost/array.hpp>
  40
  41 #include "libs/math/test/handle_test_result.hpp"
  42 #include "libs/math/test/table_type.hpp"
  43 #include "libs/math/test/functor.hpp"
  44
  45 //
  46 // Defining TEST_CPP_DEC_FLOAT enables testing of multiprecision support.
  47 // This requires the multiprecision library from sandbox/big_number.
  48 // Note that these tests *do not pass*, but they do give an idea of the
  49 // error rates that can be expected....
  50 //
  51 #ifdef TEST_CPP_DEC_FLOAT
  52 #include <boost/multiprecision/cpp_dec_float.hpp>
  53
  54 template <class R>
  55 inline R convert_to(const char* s)
  56 {
  57    try{
  58       return boost::lexical_cast<R>(s);
  59    }
  60    catch(const boost::bad_lexical_cast&)
  61    {
  62       return 0;
  63    }
  64 }
  65
  66 #define SC_(x) convert_to<T>(BOOST_STRINGIZE(x))
  67 #endif
  68
  69 #include "owens_t_T7.hpp"
  70
  71 #include <iostream>
  72 using std::cout;
  73 using std::endl;
  74 #include <limits>
  75 using std::numeric_limits;
  76
  77 void expected_results()
  78 {
  79    //
  80    // Define the max and mean errors expected for
  81    // various compilers and platforms.
  82    //
  83    const char* largest_type;
  84 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  85    if(boost::math::policies::digits<double, boost::math::policies::policy<> >() == boost::math::policies::digits<long double, boost::math::policies::policy<> >())
  86    {
  87       largest_type = "(long\\s+)?double|real_concept";
  88    }
  89    else
  90    {
  91       largest_type = "long double|real_concept";
  92    }
  93 #else
  94    largest_type = "(long\\s+)?double";
  95 #endif
  96
  97    //
  98    // Catch all cases come last:
  99    //
 100    if(std::numeric_limits<long double>::digits > 60)
 101    {
 102       add_expected_result(
 103          ".*",                            // compiler
 104          ".*",                            // stdlib
 105          ".*",                            // platform
 106          largest_type,                    // test type(s)
 107          ".*",      // test data group
 108          "boost::math::owens_t", 500, 100);  // test function
 109    }
 110    else
 111    {
 112       add_expected_result(
 113          ".*",                            // compiler
 114          ".*",                            // stdlib
 115          ".*",                            // platform
 116          largest_type,                    // test type(s)
 117          ".*",      // test data group
 118          "boost::math::owens_t", 60, 5);  // test function
 119    }
 120    //
 121    // Finish off by printing out the compiler/stdlib/platform names,
 122    // we do this to make it easier to mark up expected error rates.
 123    //
 124    std::cout << "Tests run with " << BOOST_COMPILER << ", "
 125       << BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl;
 126 }
 127
 128
 129 template <class RealType>
 130 void test_spot(
 131      RealType h,    //
 132      RealType a,    //
 133      RealType tol)   // Test tolerance
 134 {
 135   BOOST_CHECK_CLOSE_FRACTION(owens_t(h, a), 3.89119302347013668966224771378e-2L, tol);
 136 }
 137
 138
 139 template <class RealType> // Any floating-point type RealType.
 140 void test_spots(RealType)
 141 {
 142   // Basic sanity checks, test data is as accurate as long double,
 143   // so set tolerance to a few epsilon expressed as a fraction.
 144   RealType tolerance = boost::math::tools::epsilon<RealType>() * 30; // most OK with 3 eps tolerance.
 145   cout << "Tolerance = " << tolerance << "." << endl;
 146
 147   using  ::boost::math::owens_t;
 148   using ::boost::math::normal_distribution;
 149   BOOST_MATH_STD_USING // ADL of std names.
 150
 151   // Checks of six sub-methods T1 to T6.
 152   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.0625L), static_cast<RealType>(0.25L)), static_cast<RealType>(3.89119302347013668966224771378e-2L), tolerance);  // T1
 153   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(6.5L), static_cast<RealType>(0.4375L)), static_cast<RealType>(2.00057730485083154100907167685E-11L), tolerance); // T2
 154   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(7L), static_cast<RealType>( 0.96875L)), static_cast<RealType>(6.39906271938986853083219914429E-13L), tolerance); // T3
 155   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(4.78125L), static_cast<RealType>(0.0625L)), static_cast<RealType>(1.06329748046874638058307112826E-7L), tolerance); // T4
 156   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(2.L), static_cast<RealType>(0.5L)), static_cast<RealType>(8.62507798552150713113488319155E-3L), tolerance); // T5
 157   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(1.L), static_cast<RealType>(0.9999975L)), static_cast<RealType>(6.67418089782285927715589822405E-2L), tolerance); // T6
 158   //BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(L), static_cast<RealType>(L)), static_cast<RealType>(L), tolerance);
 159
 160 //   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(L), static_cast<RealType>(L)), static_cast<RealType>(L), tolerance);
 161
 162   // Spots values using Mathematica
 163   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(6.5L), static_cast<RealType>(0.4375L)), static_cast<RealType>(2.00057730485083154100907167684918851101649922551817956120806662022118024594547E-11L), tolerance);
 164   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.4375L), static_cast<RealType>(6.5L)), static_cast<RealType>(0.16540130125449396247498691826626273249659241838438244251206819782787761751256L), tolerance);
 165   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(7.L), static_cast<RealType>(0.96875L)), static_cast<RealType>(6.39906271938986853083219914428916013764797190941459233223182225724846022843930e-13L), tolerance);
 166   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.96875L), static_cast<RealType>(7.L)), static_cast<RealType>(0.08316748474602973770533230453272140919966614259525787470390475393923633179072L), tolerance);
 167   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(4.78125L), static_cast<RealType>(0.0625L)), static_cast<RealType>(1.06329748046874638058307112826015825291136503488102191050906959246644942646701e-7L), tolerance);
 168   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.0625L), static_cast<RealType>(4.78125L)), static_cast<RealType>(0.21571185819897989857261253680409017017649352928888660746045361855686569265171L), tolerance);
 169   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(2.L), static_cast<RealType>(0.5L)), static_cast<RealType>(0.00862507798552150713113488319154637187875641190390854291100809449487812876461L), tolerance);
 170   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.5L), static_cast<RealType>(2L)), static_cast<RealType>(0.14158060365397839346662819588111542648867283386549027383784843786494855594607L), tolerance);
 171
 172   // check basic properties
 173   BOOST_CHECK_EQUAL(owens_t(static_cast<RealType>(0.5L), static_cast<RealType>(2L)), owens_t(static_cast<RealType>(-0.5L), static_cast<RealType>(2L)));
 174   BOOST_CHECK_EQUAL(owens_t(static_cast<RealType>(0.5L), static_cast<RealType>(2L)), -owens_t(static_cast<RealType>(0.5L), static_cast<RealType>(-2L)));
 175   BOOST_CHECK_EQUAL(owens_t(static_cast<RealType>(0.5L), static_cast<RealType>(2L)), -owens_t(static_cast<RealType>(-0.5L), static_cast<RealType>(-2L)));
 176
 177   // Special relations from Owen's original paper:
 178   BOOST_CHECK_EQUAL(owens_t(static_cast<RealType>(0.5), static_cast<RealType>(0)), static_cast<RealType>(0));
 179   BOOST_CHECK_EQUAL(owens_t(static_cast<RealType>(10), static_cast<RealType>(0)), static_cast<RealType>(0));
 180   BOOST_CHECK_EQUAL(owens_t(static_cast<RealType>(10000), static_cast<RealType>(0)), static_cast<RealType>(0));
 181
 182   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0), static_cast<RealType>(2L)), atan(static_cast<RealType>(2L)) / (boost::math::constants::pi<RealType>() * 2), tolerance);
 183   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0), static_cast<RealType>(0.5L)), atan(static_cast<RealType>(0.5L)) / (boost::math::constants::pi<RealType>() * 2), tolerance);
 184   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0), static_cast<RealType>(2000L)), atan(static_cast<RealType>(2000L)) / (boost::math::constants::pi<RealType>() * 2), tolerance);
 185   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(5), static_cast<RealType>(1)), cdf(normal_distribution<RealType>(), 5) * cdf(complement(normal_distribution<RealType>(), 5)) / 2, tolerance);
 186   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.125), static_cast<RealType>(1)), cdf(normal_distribution<RealType>(), 0.125) * cdf(complement(normal_distribution<RealType>(), 0.125)) / 2, tolerance);
 187   if(std::numeric_limits<RealType>::has_infinity)
 188   {
 189     BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.125), std::numeric_limits<RealType>::infinity()), cdf(complement(normal_distribution<RealType>(), 0.125)) / 2, tolerance);
 190     BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(5), std::numeric_limits<RealType>::infinity()), cdf(complement(normal_distribution<RealType>(), 5)) / 2, tolerance);
 191     BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(-0.125), std::numeric_limits<RealType>::infinity()), cdf(normal_distribution<RealType>(), -0.125) / 2, tolerance);
 192     BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(-5), std::numeric_limits<RealType>::infinity()), cdf(normal_distribution<RealType>(), -5) / 2, tolerance);
 193   }
 194 } // template <class RealType>void test_spots(RealType)
 195
 196 template <class RealType> // Any floating-point type RealType.
 197 void check_against_T7(RealType)
 198 {
 199   // Basic sanity checks, test data is as accurate as long double,
 200   // so set tolerance to a few epsilon expressed as a fraction.
 201   RealType tolerance = boost::math::tools::epsilon<RealType>() * 150; // most OK with 3 eps tolerance.
 202   cout << "Tolerance = " << tolerance << "." << endl;
 203
 204   using  ::boost::math::owens_t;
 205   using namespace std; // ADL of std names.
 206
 207   // apply log scale because points near zero are more interesting
 208   for(RealType a = static_cast<RealType>(-10.0l); a < static_cast<RealType>(3l); a+= static_cast<RealType>(0.2l))
 209     for(RealType h = static_cast<RealType>(-10.0l); h < static_cast<RealType>(3.5l); h+= static_cast<RealType>(0.2l))
 210     {
 211       const RealType expa = exp(a);
 212       const RealType exph = exp(h);
 213       const RealType t = boost::math::owens_t(exph, expa);
 214       RealType t7 = boost::math::owens_t_T7(exph,expa);
 215       //if(!(boost::math::isnormal)(t) || !(boost::math::isnormal)(t7))
 216       //   std::cout << "a = " << expa << " h = " << exph << " t = " << t << " t7 = " << t7 << std::endl;
 217       BOOST_CHECK_CLOSE_FRACTION(t, t7, tolerance);
 218     }
 219
 220 } // template <class RealType>void test_spots(RealType)
 221
 222 template <class Real, class T>
 223 void do_test_owens_t(const T& data, const char* type_name, const char* test_name)
 224 {
 225    typedef typename T::value_type row_type;
 226    typedef Real                   value_type;
 227
 228    typedef value_type (*pg)(value_type, value_type);
 229 #if defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
 230    pg funcp = boost::math::owens_t<value_type>;
 231 #else
 232    pg funcp = boost::math::owens_t;
 233 #endif
 234
 235    boost::math::tools::test_result<value_type> result;
 236
 237    std::cout << "Testing " << test_name << " with type " << type_name
 238       << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
 239
 240    //
 241    // test hermite against data:
 242    //
 243    result = boost::math::tools::test_hetero<Real>(
 244       data,
 245       bind_func<Real>(funcp, 0, 1),
 246       extract_result<Real>(2));
 247    handle_test_result(result, data[result.worst()], result.worst(), type_name, "boost::math::owens_t", test_name);
 248
 249    std::cout << std::endl;
 250 }
 251
 252 template <class T>
 253 void test_owens_t(T, const char* name)
 254 {
 255    //
 256    // The actual test data is rather verbose, so it's in a separate file
 257    //
 258    // The contents are as follows, each row of data contains
 259    // three items, input value a, input value b and erf(a, b):
 260    //
 261 #  include "owens_t.ipp"
 262
 263    do_test_owens_t<T>(owens_t, name, "Owens T (medium small values)");
 264
 265 #include "owens_t_large_data.ipp"
 266
 267    do_test_owens_t<T>(owens_t_large_data, name, "Owens T (large and diverse values)");
 268 }
 269
 270
 271 BOOST_AUTO_TEST_CASE( test_main )
 272 {
 273   BOOST_MATH_CONTROL_FP;
 274
 275   expected_results();
 276
 277   // Basic sanity-check spot values.
 278
 279   // (Parameter value, arbitrarily zero, only communicates the floating point type).
 280   test_spots(0.0F); // Test float.
 281   test_spots(0.0); // Test double.
 282 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 283   test_spots(0.0L); // Test long double.
 284 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 285   test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
 286 #endif
 287 #endif
 288
 289   check_against_T7(0.0F); // Test float.
 290   check_against_T7(0.0); // Test double.
 291 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 292   check_against_T7(0.0L); // Test long double.
 293 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 294   check_against_T7(boost::math::concepts::real_concept(0.)); // Test real concept.
 295 #endif
 296 #endif
 297
 298   test_owens_t(0.0F, "float"); // Test float.
 299   test_owens_t(0.0, "double"); // Test double.
 300 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 301   test_owens_t(0.0L, "long double"); // Test long double.
 302 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 303   test_owens_t(boost::math::concepts::real_concept(0.), "real_concept"); // Test real concept.
 304 #endif
 305 #endif
 306 #ifdef TEST_CPP_DEC_FLOAT
 307   typedef boost::multiprecision::mp_number<boost::multiprecision::cpp_dec_float<35> > cpp_dec_float_35;
 308   test_owens_t(cpp_dec_float_35(0), "cpp_dec_float_35"); // Test real concept.
 309   test_owens_t(boost::multiprecision::cpp_dec_float_50(0), "cpp_dec_float_50"); // Test real concept.
 310   test_owens_t(boost::multiprecision::cpp_dec_float_100(0), "cpp_dec_float_100"); // Test real concept.
 311 #endif
 312
 313 } // BOOST_AUTO_TEST_CASE( test_main )
 314
 315 /*
 316
 317 Output:
 318
 319   Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Debug\test_owens_t.exe"
 320   Running 1 test case...
 321   Tests run with Microsoft Visual C++ version 10.0, Dinkumware standard library version 520, Win32
 322   Tolerance = 3.57628e-006.
 323   Tolerance = 6.66134e-015.
 324   Tolerance = 6.66134e-015.
 325   Tolerance = 6.66134e-015.
 326   Tolerance = 1.78814e-005.
 327   Tolerance = 3.33067e-014.
 328   Tolerance = 3.33067e-014.
 329   Tolerance = 3.33067e-014.
 330   Testing Owens T (medium small values) with type float
 331   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 332   boost::math::owens_t<float> Max = 0 RMS Mean=0
 333
 334
 335   Testing Owens T (large and diverse values) with type float
 336   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 337   boost::math::owens_t<float> Max = 0 RMS Mean=0
 338
 339
 340   Testing Owens T (medium small values) with type double
 341   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 342   boost::math::owens_t<double> Max = 4.375 RMS Mean=0.9728
 343       worst case at row: 81
 344       { 4.4206809997558594, 0.1269868016242981, 1.0900281236140834e-006 }
 345
 346
 347   Testing Owens T (large and diverse values) with type double
 348   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 349   boost::math::owens_t<double> Max = 3.781 RMS Mean=0.6206
 350       worst case at row: 430
 351       { 3.4516773223876953, 0.0010718167759478092, 4.413983645332431e-007 }
 352
 353
 354   Testing Owens T (medium small values) with type long double
 355   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 356   boost::math::owens_t<long double> Max = 4.375 RMS Mean=0.9728
 357       worst case at row: 81
 358       { 4.4206809997558594, 0.1269868016242981, 1.0900281236140834e-006 }
 359
 360
 361   Testing Owens T (large and diverse values) with type long double
 362   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 363   boost::math::owens_t<long double> Max = 3.781 RMS Mean=0.6206
 364       worst case at row: 430
 365       { 3.4516773223876953, 0.0010718167759478092, 4.413983645332431e-007 }
 366
 367
 368   Testing Owens T (medium small values) with type real_concept
 369   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 370   boost::math::owens_t<real_concept> Max = 4.375 RMS Mean=1.032
 371       worst case at row: 81
 372       { 4.4206809997558594, 0.1269868016242981, 1.0900281236140834e-006 }
 373
 374
 375   Testing Owens T (large and diverse values) with type real_concept
 376   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 377   boost::math::owens_t<real_concept> Max = 21.04 RMS Mean=1.102
 378       worst case at row: 439
 379       { 3.4516773223876953, 0.98384737968444824, 0.00013923002576038691 }
 380
 381
 382
 383   *** No errors detected
 384
 385
 386 */
 387
 388
 389