libs/math/test/test_owens_t.hpp

   1 //  (C) Copyright John Maddock 2007.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
   7 #include <boost/math/concepts/real_concept.hpp>
   8 #define BOOST_TEST_MAIN
   9 #include <boost/test/unit_test.hpp>
  10 #include <boost/test/tools/floating_point_comparison.hpp>
  11 #include <boost/math/special_functions/math_fwd.hpp>
  12 #include <boost/math/distributions/normal.hpp>
  13 #include <boost/type_traits/is_floating_point.hpp>
  14 #include <boost/array.hpp>
  15 #include "functor.hpp"
  16
  17 #include "handle_test_result.hpp"
  18 #include "table_type.hpp"
  19 #include "owens_t_T7.hpp"
  20
  21
  22 template <class RealType>
  23 void test_spot(
  24    RealType h,    //
  25    RealType a,    //
  26    RealType tol)   // Test tolerance
  27 {
  28    BOOST_CHECK_CLOSE_FRACTION(owens_t(h, a), 3.89119302347013668966224771378e-2L, tol);
  29 }
  30
  31 template <class RealType> // Any floating-point type RealType.
  32 void test_spots(RealType)
  33 {
  34    using namespace std;
  35    // Basic sanity checks, test data is as accurate as long double,
  36    // so set tolerance to a few epsilon expressed as a fraction.
  37    RealType tolerance = boost::math::tools::epsilon<RealType>() * 30; // most OK with 3 eps tolerance.
  38    cout << "Tolerance = " << tolerance << "." << endl;
  39
  40    using  ::boost::math::owens_t;
  41    using ::boost::math::normal_distribution;
  42    BOOST_MATH_STD_USING // ADL of std names.
  43
  44       // Checks of six sub-methods T1 to T6.
  45       BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.0625L), static_cast<RealType>(0.25L)), static_cast<RealType>(3.89119302347013668966224771378e-2L), tolerance);  // T1
  46    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(6.5L), static_cast<RealType>(0.4375L)), static_cast<RealType>(2.00057730485083154100907167685E-11L), tolerance); // T2
  47    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(7L), static_cast<RealType>(0.96875L)), static_cast<RealType>(6.39906271938986853083219914429E-13L), tolerance); // T3
  48    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(4.78125L), static_cast<RealType>(0.0625L)), static_cast<RealType>(1.06329748046874638058307112826E-7L), tolerance); // T4
  49    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(2.L), static_cast<RealType>(0.5L)), static_cast<RealType>(8.62507798552150713113488319155E-3L), tolerance); // T5
  50    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(1.L), static_cast<RealType>(0.9999975L)), static_cast<RealType>(6.67418089782285927715589822405E-2L), tolerance); // T6
  51    //BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(L), static_cast<RealType>(L)), static_cast<RealType>(L), tolerance);
  52
  53    //   BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(L), static_cast<RealType>(L)), static_cast<RealType>(L), tolerance);
  54
  55    // Spots values using Mathematica
  56    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(6.5L), static_cast<RealType>(0.4375L)), static_cast<RealType>(2.00057730485083154100907167684918851101649922551817956120806662022118024594547E-11L), tolerance);
  57    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.4375L), static_cast<RealType>(6.5L)), static_cast<RealType>(0.16540130125449396247498691826626273249659241838438244251206819782787761751256L), tolerance);
  58    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(7.L), static_cast<RealType>(0.96875L)), static_cast<RealType>(6.39906271938986853083219914428916013764797190941459233223182225724846022843930e-13L), tolerance);
  59    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.96875L), static_cast<RealType>(7.L)), static_cast<RealType>(0.08316748474602973770533230453272140919966614259525787470390475393923633179072L), tolerance);
  60    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(4.78125L), static_cast<RealType>(0.0625L)), static_cast<RealType>(1.06329748046874638058307112826015825291136503488102191050906959246644942646701e-7L), tolerance);
  61    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.0625L), static_cast<RealType>(4.78125L)), static_cast<RealType>(0.21571185819897989857261253680409017017649352928888660746045361855686569265171L), tolerance);
  62    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(2.L), static_cast<RealType>(0.5L)), static_cast<RealType>(0.00862507798552150713113488319154637187875641190390854291100809449487812876461L), tolerance);
  63    BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.5L), static_cast<RealType>(2L)), static_cast<RealType>(0.14158060365397839346662819588111542648867283386549027383784843786494855594607L), tolerance);
  64
  65    // check basic properties
  66    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)));
  67    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)));
  68    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)));
  69
  70    // Special relations from Owen's original paper:
  71    BOOST_CHECK_EQUAL(owens_t(static_cast<RealType>(0.5), static_cast<RealType>(0)), static_cast<RealType>(0));
  72    BOOST_CHECK_EQUAL(owens_t(static_cast<RealType>(10), static_cast<RealType>(0)), static_cast<RealType>(0));
  73    BOOST_CHECK_EQUAL(owens_t(static_cast<RealType>(10000), static_cast<RealType>(0)), static_cast<RealType>(0));
  74
  75    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);
  76    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);
  77    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);
  78    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);
  79    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);
  80    if(std::numeric_limits<RealType>::has_infinity)
  81    {
  82       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);
  83       BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(5), std::numeric_limits<RealType>::infinity()), cdf(complement(normal_distribution<RealType>(), 5)) / 2, tolerance);
  84       BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(-0.125), std::numeric_limits<RealType>::infinity()), cdf(normal_distribution<RealType>(), -0.125) / 2, tolerance);
  85       BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(-5), std::numeric_limits<RealType>::infinity()), cdf(normal_distribution<RealType>(), -5) / 2, tolerance);
  86    }
  87 } // template <class RealType>void test_spots(RealType)
  88
  89 template <class RealType> // Any floating-point type RealType.
  90 void check_against_T7(RealType)
  91 {
  92    using namespace std;
  93    // Basic sanity checks, test data is as accurate as long double,
  94    // so set tolerance to a few epsilon expressed as a fraction.
  95    RealType tolerance = boost::math::tools::epsilon<RealType>() * 150; // most OK with 3 eps tolerance.
  96    cout << "Tolerance = " << tolerance << "." << endl;
  97
  98    using  ::boost::math::owens_t;
  99    using namespace std; // ADL of std names.
 100
 101    // apply log scale because points near zero are more interesting
 102    for(RealType a = static_cast<RealType>(-10.0l); a < static_cast<RealType>(3l); a += static_cast<RealType>(0.2l))
 103       for(RealType h = static_cast<RealType>(-10.0l); h < static_cast<RealType>(3.5l); h += static_cast<RealType>(0.2l))
 104       {
 105          const RealType expa = exp(a);
 106          const RealType exph = exp(h);
 107          const RealType t = boost::math::owens_t(exph, expa);
 108          RealType t7 = boost::math::owens_t_T7(exph, expa);
 109          //if(!(boost::math::isnormal)(t) || !(boost::math::isnormal)(t7))
 110          //   std::cout << "a = " << expa << " h = " << exph << " t = " << t << " t7 = " << t7 << std::endl;
 111          BOOST_CHECK_CLOSE_FRACTION(t, t7, tolerance);
 112       }
 113
 114 } // template <class RealType>void test_spots(RealType)
 115
 116 template <class Real, class T>
 117 void do_test_owens_t(const T& data, const char* type_name, const char* test_name)
 118 {
 119 #if !(defined(ERROR_REPORTING_MODE) && !defined(OWENS_T_FUNCTION_TO_TEST))
 120    typedef Real                   value_type;
 121
 122    typedef value_type(*pg)(value_type, value_type);
 123 #ifdef OWENS_T_FUNCTION_TO_TEST
 124    pg funcp = OWENS_T_FUNCTION_TO_TEST;
 125 #elif defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS)
 126    pg funcp = boost::math::owens_t<value_type>;
 127 #else
 128    pg funcp = boost::math::owens_t;
 129 #endif
 130
 131    boost::math::tools::test_result<value_type> result;
 132
 133    std::cout << "Testing " << test_name << " with type " << type_name
 134       << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";
 135
 136    //
 137    // test owens_t against data:
 138    //
 139    result = boost::math::tools::test_hetero<Real>(
 140       data,
 141       bind_func<Real>(funcp, 0, 1),
 142       extract_result<Real>(2));
 143    handle_test_result(result, data[result.worst()], result.worst(), type_name, "owens_t", test_name);
 144
 145    std::cout << std::endl;
 146 #endif
 147 }
 148
 149 template <class T>
 150 void test_owens_t(T, const char* name)
 151 {
 152    //
 153    // The actual test data is rather verbose, so it's in a separate file
 154    //
 155    // The contents are as follows, each row of data contains
 156    // three items, input value a, input value b and erf(a, b):
 157    //
 158 #  include "owens_t.ipp"
 159
 160    do_test_owens_t<T>(owens_t, name, "Owens T (medium small values)");
 161
 162 #include "owens_t_large_data.ipp"
 163
 164    do_test_owens_t<T>(owens_t_large_data, name, "Owens T (large and diverse values)");
 165 }