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