libs/math/test/test_nc_beta.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 #ifndef BOOST_MATH_OVERFLOW_ERROR_POLICY
   7 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
   8 #endif
   9 #include <boost/math/concepts/real_concept.hpp>
  10 #define BOOST_TEST_MAIN
  11 #include <boost/test/unit_test.hpp>
  12 #include <boost/test/tools/floating_point_comparison.hpp>
  13 #include <boost/math/distributions/non_central_beta.hpp>
  14 #include <boost/math/distributions/poisson.hpp>
  15 #include <boost/type_traits/is_floating_point.hpp>
  16 #include <boost/array.hpp>
  17 #include "functor.hpp"
  18
  19 #include "handle_test_result.hpp"
  20 #include "table_type.hpp"
  21
  22 #define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \
  23       {\
  24       unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
  25       BOOST_CHECK_CLOSE(a, b, prec); \
  26       if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
  27             {\
  28          std::cerr << "Failure was at row " << i << std::endl;\
  29          std::cerr << std::setprecision(35); \
  30          std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
  31          std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
  32             }\
  33       }
  34
  35 #define BOOST_CHECK_EX(a, i) \
  36       {\
  37       unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
  38       BOOST_CHECK(a); \
  39       if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
  40             {\
  41          std::cerr << "Failure was at row " << i << std::endl;\
  42          std::cerr << std::setprecision(35); \
  43          std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
  44          std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
  45             }\
  46       }
  47
  48 template <class T>
  49 T nc_beta_cdf(T a, T b, T nc, T x)
  50 {
  51 #ifdef NC_BETA_CDF_FUNCTION_TO_TEST
  52    return NC_BETA_CDF_FUNCTION_TO_TEST(a, b, nc, x);
  53 #else
  54    return cdf(boost::math::non_central_beta_distribution<T>(a, b, nc), x);
  55 #endif
  56 }
  57
  58 template <class T>
  59 T nc_beta_ccdf(T a, T b, T nc, T x)
  60 {
  61 #ifdef NC_BETA_CCDF_FUNCTION_TO_TEST
  62    return NC_BETA_CCDF_FUNCTION_TO_TEST(a, b, nc, x);
  63 #else
  64    return cdf(complement(boost::math::non_central_beta_distribution<T>(a, b, nc), x));
  65 #endif
  66 }
  67
  68 template <typename Real, typename T>
  69 void do_test_nc_chi_squared(T& data, const char* type_name, const char* test)
  70 {
  71    typedef Real                   value_type;
  72
  73    std::cout << "Testing: " << test << std::endl;
  74
  75    value_type(*fp1)(value_type, value_type, value_type, value_type) = nc_beta_cdf;
  76    boost::math::tools::test_result<value_type> result;
  77
  78 #if !(defined(ERROR_REPORTING_MODE) && !defined(NC_BETA_CDF_FUNCTION_TO_TEST))
  79    result = boost::math::tools::test_hetero<Real>(
  80       data,
  81       bind_func<Real>(fp1, 0, 1, 2, 3),
  82       extract_result<Real>(4));
  83    handle_test_result(result, data[result.worst()], result.worst(),
  84       type_name, "non central beta CDF", test);
  85 #endif
  86 #if !(defined(ERROR_REPORTING_MODE) && !defined(NC_BETA_CCDF_FUNCTION_TO_TEST))
  87    fp1 = nc_beta_ccdf;
  88    result = boost::math::tools::test_hetero<Real>(
  89       data,
  90       bind_func<Real>(fp1, 0, 1, 2, 3),
  91       extract_result<Real>(5));
  92    handle_test_result(result, data[result.worst()], result.worst(),
  93       type_name, "non central beta CDF complement", test);
  94 #endif
  95    std::cout << std::endl;
  96
  97 }
  98
  99 template <typename Real, typename T>
 100 void quantile_sanity_check(T& data, const char* type_name, const char* test)
 101 {
 102 #ifndef ERROR_REPORTING_MODE
 103    typedef Real                   value_type;
 104
 105    //
 106    // Tests with type real_concept take rather too long to run, so
 107    // for now we'll disable them:
 108    //
 109    if(!boost::is_floating_point<value_type>::value)
 110       return;
 111
 112    std::cout << "Testing: " << type_name << " quantile sanity check, with tests " << test << std::endl;
 113
 114    //
 115    // These sanity checks test for a round trip accuracy of one half
 116    // of the bits in T, unless T is type float, in which case we check
 117    // for just one decimal digit.  The problem here is the sensitivity
 118    // of the functions, not their accuracy.  This test data was generated
 119    // for the forward functions, which means that when it is used as
 120    // the input to the inverses then it is necessarily inexact.  This rounding
 121    // of the input is what makes the data unsuitable for use as an accuracy check,
 122    // and also demonstrates that you can't in general round-trip these functions.
 123    // It is however a useful sanity check.
 124    //
 125    value_type precision = static_cast<value_type>(ldexp(1.0, 1 - boost::math::policies::digits<value_type, boost::math::policies::policy<> >() / 2)) * 100;
 126    if(boost::math::policies::digits<value_type, boost::math::policies::policy<> >() < 50)
 127       precision = 1;   // 1% or two decimal digits, all we can hope for when the input is truncated to float
 128
 129    for(unsigned i = 0; i < data.size(); ++i)
 130    {
 131       //
 132       // Test case 493 fails at float precision: not enough bits to get
 133       // us back where we started:
 134       //
 135       if((i == 493) && boost::is_same<float, value_type>::value)
 136          continue;
 137
 138       if(data[i][4] == 0)
 139       {
 140          BOOST_CHECK(0 == quantile(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1], data[i][2]), data[i][4]));
 141       }
 142       else if(data[i][4] < 0.9999f)
 143       {
 144          value_type p = quantile(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1], data[i][2]), data[i][4]);
 145          value_type pt = data[i][3];
 146          BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
 147       }
 148       if(data[i][5] == 0)
 149       {
 150          BOOST_CHECK(1 == quantile(complement(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1], data[i][2]), data[i][5])));
 151       }
 152       else if(data[i][5] < 0.9999f)
 153       {
 154          value_type p = quantile(complement(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1], data[i][2]), data[i][5]));
 155          value_type pt = data[i][3];
 156          BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
 157       }
 158       if(boost::math::tools::digits<value_type>() > 50)
 159       {
 160          //
 161          // Sanity check mode, accuracy of
 162          // the mode is at *best* the square root of the accuracy of the PDF:
 163          //
 164          value_type m = mode(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1], data[i][2]));
 165          if((m == 1) || (m == 0))
 166             break;
 167          value_type p = pdf(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1], data[i][2]), m);
 168          if(m * (1 + sqrt(precision) * 10) < 1)
 169          {
 170             BOOST_CHECK_EX(pdf(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1], data[i][2]), m * (1 + sqrt(precision) * 10)) <= p, i);
 171          }
 172          if(m * (1 - sqrt(precision)) * 10 > boost::math::tools::min_value<value_type>())
 173          {
 174             BOOST_CHECK_EX(pdf(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1], data[i][2]), m * (1 - sqrt(precision)) * 10) <= p, i);
 175          }
 176       }
 177    }
 178 #endif
 179 }
 180
 181 template <typename T>
 182 void test_accuracy(T, const char* type_name)
 183 {
 184 #if !defined(TEST_DATA) || (TEST_DATA == 1)
 185 #include "ncbeta.ipp"
 186    do_test_nc_chi_squared<T>(ncbeta, type_name, "Non Central Beta, medium parameters");
 187    quantile_sanity_check<T>(ncbeta, type_name, "Non Central Beta, medium parameters");
 188 #endif
 189 #if !defined(TEST_DATA) || (TEST_DATA == 2)
 190 #include "ncbeta_big.ipp"
 191    do_test_nc_chi_squared<T>(ncbeta_big, type_name, "Non Central Beta, large parameters");
 192    // Takes too long to run:
 193    // quantile_sanity_check(ncbeta_big, type_name, "Non Central Beta, large parameters");
 194 #endif
 195 }