libs/math/test/test_nc_beta.cpp

   1 // test_nc_beta.cpp
   2
   3 // Copyright John Maddock 2008.
   4
   5 // Use, modification and distribution are subject to the
   6 // Boost Software License, Version 1.0.
   7 // (See accompanying file LICENSE_1_0.txt
   8 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   9
  10 //
  11 // This must appear *before* any #includes, and precludes pch usage:
  12 //
  13 #define BOOST_MATH_ASSERT_UNDEFINED_POLICY false
  14
  15 #ifdef _MSC_VER
  16 #pragma warning (disable:4127 4512)
  17 #endif
  18
  19 #if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
  20 #  define TEST_FLOAT
  21 #  define TEST_DOUBLE
  22 #  define TEST_LDOUBLE
  23 #  define TEST_REAL_CONCEPT
  24 #endif
  25
  26 #include <boost/math/concepts/real_concept.hpp> // for real_concept
  27 #include <boost/math/distributions/non_central_beta.hpp> // for chi_squared_distribution
  28 #include <boost/math/distributions/poisson.hpp> // for poisson_distribution
  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 "test_ncbeta_hooks.hpp"
  38 #include "table_type.hpp"
  39 #include "test_nc_beta.hpp"
  40
  41 #include <iostream>
  42 using std::cout;
  43 using std::endl;
  44 #include <limits>
  45 using std::numeric_limits;
  46
  47 void expected_results()
  48 {
  49    //
  50    // Define the max and mean errors expected for
  51    // various compilers and platforms.
  52    //
  53    const char* largest_type;
  54 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  55    if(boost::math::policies::digits<double, boost::math::policies::policy<> >() == boost::math::policies::digits<long double, boost::math::policies::policy<> >())
  56    {
  57       largest_type = "(long\\s+)?double|real_concept";
  58    }
  59    else
  60    {
  61       largest_type = "long double|real_concept";
  62    }
  63 #else
  64    largest_type = "(long\\s+)?double|real_concept";
  65 #endif
  66
  67 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  68    if(boost::math::tools::digits<long double>() == 64)
  69    {
  70       //
  71       // Allow a small amount of error leakage from long double to double:
  72       //
  73       add_expected_result(
  74          "[^|]*",                          // compiler
  75          "[^|]*",                          // stdlib
  76          "[^|]*",                          // platform
  77          "double",                         // test type(s)
  78          "[^|]*large[^|]*",                // test data group
  79          "[^|]*", 5, 5);                   // test function
  80    }
  81
  82    if(boost::math::tools::digits<long double>() == 64)
  83    {
  84       add_expected_result(
  85          "[^|]*",                          // compiler
  86          "[^|]*",                          // stdlib
  87          "[^|]*",                          // platform
  88          largest_type,                     // test type(s)
  89          "[^|]*medium[^|]*",               // test data group
  90          "[^|]*", 1200, 500);               // test function
  91       add_expected_result(
  92          "[^|]*",                          // compiler
  93          "[^|]*",                          // stdlib
  94          "[^|]*",                          // platform
  95          largest_type,                     // test type(s)
  96          "[^|]*large[^|]*",                // test data group
  97          "[^|]*", 40000, 6000);            // test function
  98    }
  99 #endif
 100    //
 101    // Catch all cases come last:
 102    //
 103    add_expected_result(
 104       "[^|]*",                          // compiler
 105       "[^|]*",                          // stdlib
 106       "[^|]*",                          // platform
 107       largest_type,                     // test type(s)
 108       "[^|]*medium[^|]*",               // test data group
 109       "[^|]*", 1500, 500);               // test function
 110    add_expected_result(
 111       "[^|]*",                          // compiler
 112       "[^|]*",                          // stdlib
 113       "[^|]*",                          // platform
 114       "real_concept",                   // test type(s)
 115       "[^|]*large[^|]*",                // test data group
 116       "[^|]*", 30000, 4000);             // test function
 117    add_expected_result(
 118       "[^|]*",                          // compiler
 119       "[^|]*",                          // stdlib
 120       "[^|]*",                          // platform
 121       largest_type,                     // test type(s)
 122       "[^|]*large[^|]*",                // test data group
 123       "[^|]*", 20000, 2000);             // test function
 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 template <class RealType>
 133 RealType naive_pdf(RealType a, RealType b, RealType lam, RealType x)
 134 {
 135    using namespace boost::math;
 136
 137    RealType term = pdf(poisson_distribution<RealType>(lam/2), 0)
 138       * ibeta_derivative(a, b, x);
 139    RealType sum = term;
 140
 141    int i = 1;
 142    while(term / sum > tools::epsilon<RealType>())
 143    {
 144       term = pdf(poisson_distribution<RealType>(lam/2), i)
 145       * ibeta_derivative(a + i, b, x);
 146       ++i;
 147       sum += term;
 148    }
 149    return sum;
 150 }
 151
 152 template <class RealType>
 153 void test_spot(
 154      RealType a,     // alpha
 155      RealType b,     // beta
 156      RealType ncp,   // non-centrality param
 157      RealType cs,    // Chi Square statistic
 158      RealType P,     // CDF
 159      RealType Q,     // Complement of CDF
 160      RealType D,     // PDF
 161      RealType tol)   // Test tolerance
 162 {
 163    boost::math::non_central_beta_distribution<RealType> dist(a, b, ncp);
 164    BOOST_CHECK_CLOSE(
 165       cdf(dist, cs), P, tol);
 166    //
 167    // Sanity checking using the naive PDF calculation above fails at
 168    // float precision:
 169    //
 170    if(!boost::is_same<float, RealType>::value)
 171    {
 172       BOOST_CHECK_CLOSE(
 173          pdf(dist, cs), naive_pdf(dist.alpha(), dist.beta(), ncp, cs), tol);
 174    }
 175    BOOST_CHECK_CLOSE(
 176       pdf(dist, cs), D, tol);
 177
 178    if((P < 0.99) && (Q < 0.99))
 179    {
 180       //
 181       // We can only check this if P is not too close to 1,
 182       // so that we can guarantee Q is reasonably free of error:
 183       //
 184       BOOST_CHECK_CLOSE(
 185          cdf(complement(dist, cs)), Q, tol);
 186       BOOST_CHECK_CLOSE(
 187             quantile(dist, P), cs, tol * 10);
 188       BOOST_CHECK_CLOSE(
 189             quantile(complement(dist, Q)), cs, tol * 10);
 190    }
 191 }
 192
 193 template <class RealType> // Any floating-point type RealType.
 194 void test_spots(RealType)
 195 {
 196    RealType tolerance = (std::max)(
 197       boost::math::tools::epsilon<RealType>() * 100,
 198       (RealType)1e-6) * 100;
 199    RealType abs_tolerance = boost::math::tools::epsilon<RealType>() * 100;
 200
 201    cout << "Tolerance = " << tolerance << "%." << endl;
 202
 203    //
 204    // Spot tests use values computed by the R statistical
 205    // package and the pbeta and dbeta functions:
 206    //
 207    test_spot(
 208      RealType(2),                   // alpha
 209      RealType(5),                   // beta
 210      RealType(1),                   // non-centrality param
 211      RealType(0.25),                // Chi Square statistic
 212      RealType(0.3658349),           // CDF
 213      RealType(1-0.3658349),         // Complement of CDF
 214      RealType(2.184465),            // PDF
 215      RealType(tolerance));
 216    test_spot(
 217      RealType(20),                  // alpha
 218      RealType(15),                  // beta
 219      RealType(35),                  // non-centrality param
 220      RealType(0.75),                // Chi Square statistic
 221      RealType(0.6994175),           // CDF
 222      RealType(1-0.6994175),         // Complement of CDF
 223      RealType(5.576146),            // PDF
 224      RealType(tolerance));
 225    test_spot(
 226      RealType(100),                 // alpha
 227      RealType(3),                   // beta
 228      RealType(63),                  // non-centrality param
 229      RealType(0.95),                // Chi Square statistic
 230      RealType(0.03529306),          // CDF
 231      RealType(1-0.03529306),        // Complement of CDF
 232      RealType(3.637894),            // PDF
 233      RealType(tolerance));
 234    test_spot(
 235      RealType(0.25),                // alpha
 236      RealType(0.75),                // beta
 237      RealType(150),                 // non-centrality param
 238      RealType(0.975),               // Chi Square statistic
 239      RealType(0.09752216),          // CDF
 240      RealType(1-0.09752216),        // Complement of CDF
 241      RealType(8.020935),            // PDF
 242      RealType(tolerance));
 243
 244    BOOST_MATH_STD_USING
 245    boost::math::non_central_beta_distribution<RealType> dist(100, 3, 63);
 246    BOOST_CHECK_CLOSE(mean(dist), RealType(4.82280451915522329944315287538684030781836554279474240490936e13L) * exp(-RealType(31.5)) * 100 / 103, tolerance);
 247    // Variance only guarantees small absolute error:
 248    BOOST_CHECK_SMALL(variance(dist)
 249       - static_cast<RealType>(RealType(4.85592267707818899235900237275021938334418424134218087127572e13L)
 250       * exp(RealType(-31.5)) * 100 * 101 / (103 * 104) -
 251       RealType(4.82280451915522329944315287538684030781836554279474240490936e13L) * RealType(4.82280451915522329944315287538684030781836554279474240490936e13L)
 252       * exp(RealType(-63)) * 10000 / (103 * 103)), abs_tolerance);
 253    BOOST_MATH_CHECK_THROW(skewness(dist), boost::math::evaluation_error);
 254    BOOST_MATH_CHECK_THROW(kurtosis(dist), boost::math::evaluation_error);
 255    BOOST_MATH_CHECK_THROW(kurtosis_excess(dist), boost::math::evaluation_error);
 256 } // template <class RealType>void test_spots(RealType)
 257
 258
 259 BOOST_AUTO_TEST_CASE( test_main )
 260 {
 261    BOOST_MATH_CONTROL_FP;
 262    // Basic sanity-check spot values.
 263     expected_results();
 264    // (Parameter value, arbitrarily zero, only communicates the floating point type).
 265 #ifdef TEST_FLOAT
 266    test_spots(0.0F); // Test float.
 267 #endif
 268 #ifdef TEST_DOUBLE
 269    test_spots(0.0); // Test double.
 270 #endif
 271 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 272 #ifdef TEST_LDOUBLE
 273    test_spots(0.0L); // Test long double.
 274 #endif
 275 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 276 #ifdef TEST_REAL_CONCEPT
 277    test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
 278 #endif
 279 #endif
 280 #endif
 281
 282 #ifdef TEST_FLOAT
 283    test_accuracy(0.0F, "float"); // Test float.
 284 #endif
 285 #ifdef TEST_DOUBLE
 286    test_accuracy(0.0, "double"); // Test double.
 287 #endif
 288 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 289 #ifdef TEST_LDOUBLE
 290    test_accuracy(0.0L, "long double"); // Test long double.
 291 #endif
 292 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 293 #ifdef TEST_REAL_CONCEPT
 294    test_accuracy(boost::math::concepts::real_concept(0.), "real_concept"); // Test real concept.
 295 #endif
 296 #endif
 297 #endif
 298
 299 } // BOOST_AUTO_TEST_CASE( test_main )
 300