libs/math/test/test_nc_f.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 #include <pch.hpp>
  11
  12 #ifdef _MSC_VER
  13 #pragma warning (disable:4127 4512)
  14 #endif
  15
  16 #if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
  17 #  define TEST_FLOAT
  18 #  define TEST_DOUBLE
  19 #  define TEST_LDOUBLE
  20 #  define TEST_REAL_CONCEPT
  21 #endif
  22
  23 #include <boost/math/concepts/real_concept.hpp> // for real_concept
  24 #include <boost/math/distributions/non_central_f.hpp> // for chi_squared_distribution
  25 #define BOOST_TEST_MAIN
  26 #include <boost/test/unit_test.hpp> // for test_main
  27 #include <boost/test/results_collector.hpp>
  28 #include <boost/test/unit_test.hpp>
  29 #include <boost/test/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE
  30 #include "test_out_of_range.hpp"
  31
  32 #include "functor.hpp"
  33 #include "handle_test_result.hpp"
  34
  35 #include <iostream>
  36 #include <iomanip>
  37 using std::cout;
  38 using std::endl;
  39 #include <limits>
  40 using std::numeric_limits;
  41
  42 #define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \
  43    {\
  44       unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
  45       BOOST_CHECK_CLOSE(a, b, prec); \
  46       if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
  47       {\
  48          std::cerr << "Failure was at row " << i << std::endl;\
  49          std::cerr << std::setprecision(35); \
  50          std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
  51          std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
  52       }\
  53    }
  54
  55 #define BOOST_CHECK_EX(a, i) \
  56    {\
  57       unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
  58       BOOST_CHECK(a); \
  59       if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
  60       {\
  61          std::cerr << "Failure was at row " << i << std::endl;\
  62          std::cerr << std::setprecision(35); \
  63          std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
  64          std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
  65       }\
  66    }
  67
  68 void expected_results()
  69 {
  70    //
  71    // Define the max and mean errors expected for
  72    // various compilers and platforms.
  73    //
  74    const char* largest_type;
  75 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  76    if(boost::math::policies::digits<double, boost::math::policies::policy<> >() == boost::math::policies::digits<long double, boost::math::policies::policy<> >())
  77    {
  78       largest_type = "(long\\s+)?double|real_concept";
  79    }
  80    else
  81    {
  82       largest_type = "long double|real_concept";
  83    }
  84 #else
  85    largest_type = "(long\\s+)?double|real_concept";
  86 #endif
  87
  88    //
  89    // Finish off by printing out the compiler/stdlib/platform names,
  90    // we do this to make it easier to mark up expected error rates.
  91    //
  92    std::cout << "Tests run with " << BOOST_COMPILER << ", "
  93       << BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl;
  94 }
  95
  96 template <class RealType>
  97 RealType naive_pdf(RealType v1, RealType v2, RealType l, RealType f)
  98 {
  99    BOOST_MATH_STD_USING
 100    RealType sum = 0;
 101    for(int i = 0; ; ++i)
 102    {
 103       RealType term = -l/2 + log(l/2) * i;
 104       term += boost::math::lgamma(v2/2 + v1/2+i) - (boost::math::lgamma(v2/2) + boost::math::lgamma(v1/2+i));
 105       term -= boost::math::lgamma(RealType(i+1));
 106       term += log(v1/v2) * (v1/2+i) + log(v2 / (v2 + v1 * f)) * ((v1 + v2) / 2 + i);
 107       term += log(f) * (v1/2 - 1 + i);
 108       term = exp(term);
 109       sum += term;
 110       if((term/sum < boost::math::tools::epsilon<RealType>()) || (term == 0))
 111          break;
 112    }
 113    return sum;
 114 }
 115
 116 template <class RealType>
 117 void test_spot(
 118      RealType a,     // df1
 119      RealType b,     // df2
 120      RealType ncp,   // non-centrality param
 121      RealType x,     // F statistic
 122      RealType P,     // CDF
 123      RealType Q,     // Complement of CDF
 124      RealType D,     // PDF
 125      RealType tol)   // Test tolerance
 126 {
 127    boost::math::non_central_f_distribution<RealType> dist(a, b, ncp);
 128    BOOST_CHECK_CLOSE(
 129       cdf(dist, x), P, tol);
 130    BOOST_CHECK_CLOSE(
 131       pdf(dist, x), D, tol);
 132    if(boost::math::tools::digits<RealType>() > 50)
 133    {
 134       //
 135       // The "naive" pdf calculation fails at float precision.
 136       //
 137       BOOST_CHECK_CLOSE(
 138          pdf(dist, x), naive_pdf(a, b, ncp, x), tol);
 139    }
 140
 141    if((P < 0.99) && (Q < 0.99))
 142    {
 143       //
 144       // We can only check this if P is not too close to 1,
 145       // so that we can guarentee Q is reasonably free of error:
 146       //
 147       BOOST_CHECK_CLOSE(
 148          cdf(complement(dist, x)), Q, tol);
 149       BOOST_CHECK_CLOSE(
 150             quantile(dist, P), x, tol * 10);
 151       BOOST_CHECK_CLOSE(
 152             quantile(complement(dist, Q)), x, tol * 10);
 153    }
 154    if(boost::math::tools::digits<RealType>() > 50)
 155    {
 156       //
 157       // Sanity check mode:
 158       //
 159       RealType m = mode(dist);
 160       RealType p = pdf(dist, m);
 161       BOOST_CHECK(pdf(dist, m * (1 + sqrt(tol) * 10)) <= p);
 162       BOOST_CHECK(pdf(dist, m * (1 - sqrt(tol)) * 10) <= p);
 163    }
 164 }
 165
 166 template <class RealType> // Any floating-point type RealType.
 167 void test_spots(RealType)
 168 {
 169    RealType tolerance = boost::math::tools::epsilon<RealType>() * 10000;
 170
 171    cout << "Tolerance = " << (tolerance / 100) << "%." << endl;
 172
 173    //
 174    // Spot tests from Mathematica computed values:
 175    //
 176    test_spot(
 177       RealType(5),                   // alpha
 178       RealType(2),                   // beta
 179       RealType(1),                   // non-centrality param
 180       RealType(1.5),                   // F statistic
 181       RealType(0.49845842011686358665786775091245664L),           // CDF
 182       RealType(1 - 0.49845842011686358665786775091245664L),         // Complement of CDF
 183       RealType(0.20251311620629730205859816288225385L),           // PDF
 184       RealType(tolerance));
 185
 186    test_spot(
 187      RealType(2),                   // alpha
 188      RealType(5),                   // beta
 189      RealType(1),                   // non-centrality param
 190      RealType(2),                   // F statistic
 191      RealType(0.64938711196845800322066756609406894L),           // CDF
 192      RealType(1 - 0.64938711196845800322066756609406894L),         // Complement of CDF
 193      RealType(0.15512617916132011524583796078456003L),           // PDF
 194      RealType(tolerance));
 195    test_spot(
 196      RealType(100),                 // alpha
 197      RealType(5),                   // beta
 198      RealType(15),                  // non-centrality param
 199      RealType(105),                 // F statistic
 200      RealType(0.99996207325249555786258005958906310L),            // CDF
 201      RealType(0.000037926747504442137419940410936905407L),          // Complement of CDF
 202      RealType(8.9562292619539161551049126260104435e-7),         // PDF
 203      RealType(tolerance * 10));
 204    test_spot(
 205      RealType(100),                 // alpha
 206      RealType(5),                   // beta
 207      RealType(15),                  // non-centrality param
 208      RealType(1.5),                 // F statistic
 209      RealType(0.57592315596686179870591317303126895L),           // CDF
 210      RealType(1 - 0.57592315596686179870591317303126895L),         // Complement of CDF
 211      RealType(0.36743745541686900593212039288061162L),           // PDF
 212      RealType(tolerance * 5));
 213    test_spot(
 214      RealType(5),                   // alpha
 215      RealType(100),                 // beta
 216      RealType(102),                 // non-centrality param
 217      RealType(25),                  // F statistic
 218      RealType(0.74993383259829917593356265102363267L),           // CDF
 219      RealType(1 - 0.74993383259829917593356265102363267L),         // Complement of CDF
 220      RealType(0.054467600423154020554785779421659007L),           // PDF
 221      RealType(tolerance * 5));
 222    test_spot(
 223      RealType(85),                  // alpha
 224      RealType(100),                 // beta
 225      RealType(0.5),                 // non-centrality param
 226      RealType(1.25),                // F statistic
 227      RealType(0.85228624977948142884820398473385455L),           // CDF
 228      RealType(0.14771375022051857115179601526614545L),         // Complement of CDF
 229      RealType(0.88510283331631848299643323511414868L),           // PDF
 230      RealType(tolerance * 5));
 231
 232    //
 233    // Spot tests use values computed by the R statistical
 234    // package and the pbeta and dbeta functions:
 235    //
 236    tolerance = (std::max)(
 237       boost::math::tools::epsilon<RealType>() * 100,
 238       (RealType)1e-6) * 100;
 239
 240    test_spot(
 241       RealType(85),                  // alpha
 242       RealType(100),                 // beta
 243       RealType(245),                 // non-centrality param
 244       RealType(3.5),                 // F statistic
 245       RealType(0.2697244),           // CDF
 246       RealType(1 - 0.2697244),         // Complement of CDF
 247       RealType(0.54352369104452836465948073900030320L),           // PDF
 248       RealType(tolerance));
 249
 250    BOOST_MATH_STD_USING
 251
 252    //
 253    // 5 eps expressed as a persentage, otherwise the limit of the test data:
 254    //
 255    RealType tol2 = (std::max)(boost::math::tools::epsilon<RealType>() * 500, RealType(1e-25));
 256    RealType x = 2;
 257
 258    boost::math::non_central_f_distribution<RealType> dist(20, 15, 30);
 259    // mean:
 260    BOOST_CHECK_CLOSE(
 261       mean(dist)
 262       , static_cast<RealType>(2.8846153846153846153846153846154L), tol2);
 263    // variance:
 264    BOOST_CHECK_CLOSE(
 265       variance(dist)
 266       , static_cast<RealType>(2.1422807961269499731038192576654L), tol2);
 267    // std deviation:
 268    BOOST_CHECK_CLOSE(
 269       standard_deviation(dist)
 270       , sqrt(variance(dist)), tol2);
 271    // hazard:
 272    BOOST_CHECK_CLOSE(
 273       hazard(dist, x)
 274       , pdf(dist, x) / cdf(complement(dist, x)), tol2);
 275    // cumulative hazard:
 276    BOOST_CHECK_CLOSE(
 277       chf(dist, x)
 278       , -log(cdf(complement(dist, x))), tol2);
 279    // coefficient_of_variation:
 280    BOOST_CHECK_CLOSE(
 281       coefficient_of_variation(dist)
 282       , standard_deviation(dist) / mean(dist), tol2);
 283    BOOST_CHECK_CLOSE(
 284       median(dist),
 285       quantile(
 286       dist,
 287       static_cast<RealType>(0.5)), static_cast<RealType>(tol2));
 288    // mode:
 289    BOOST_CHECK_CLOSE(
 290       mode(dist)
 291       , static_cast<RealType>(2.070019130232759428074835788815387743293972985542L), sqrt(tolerance));
 292    // skewness:
 293    BOOST_CHECK_CLOSE(
 294       skewness(dist)
 295       , static_cast<RealType>(2.1011821125804540669752380228510691320707051692719L), tol2);
 296    // kurtosis:
 297    BOOST_CHECK_CLOSE(
 298       kurtosis(dist)
 299       , 3 + kurtosis_excess(dist), tol2);
 300    // kurtosis excess:
 301    BOOST_CHECK_CLOSE(
 302       kurtosis_excess(dist)
 303       , static_cast<RealType>(13.225781681053154767604638331440974359675882226873L), tol2);
 304
 305    // Error handling checks:
 306    check_out_of_range<boost::math::non_central_f_distribution<RealType> >(1, 1, 1);
 307    BOOST_CHECK_THROW(pdf(boost::math::non_central_f_distribution<RealType>(0, 1, 1), 0), std::domain_error);
 308    BOOST_CHECK_THROW(pdf(boost::math::non_central_f_distribution<RealType>(-1, 1, 1), 0), std::domain_error);
 309    BOOST_CHECK_THROW(pdf(boost::math::non_central_f_distribution<RealType>(1, 0, 1), 0), std::domain_error);
 310    BOOST_CHECK_THROW(pdf(boost::math::non_central_f_distribution<RealType>(1, -1, 1), 0), std::domain_error);
 311    BOOST_CHECK_THROW(quantile(boost::math::non_central_f_distribution<RealType>(1, 1, 1), -1), std::domain_error);
 312    BOOST_CHECK_THROW(quantile(boost::math::non_central_f_distribution<RealType>(1, 1, 1), 2), std::domain_error);
 313 } // template <class RealType>void test_spots(RealType)
 314
 315 BOOST_AUTO_TEST_CASE( test_main )
 316 {
 317    BOOST_MATH_CONTROL_FP;
 318    // Basic sanity-check spot values.
 319    expected_results();
 320    // (Parameter value, arbitrarily zero, only communicates the floating point type).
 321    test_spots(0.0F); // Test float.
 322    test_spots(0.0); // Test double.
 323 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 324    test_spots(0.0L); // Test long double.
 325 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 326    test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
 327 #endif
 328 #endif
 329
 330
 331 } // BOOST_AUTO_TEST_CASE( test_main )
 332