libs/math/test/test_exponential_dist.cpp

   1 // Copyright John Maddock 2006.
   2 // Copyright Paul A. Bristow 2007.
   3
   4 // Use, modification and distribution are subject to the
   5 // Boost Software License, Version 1.0.
   6 // (See accompanying file LICENSE_1_0.txt
   7 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   8
   9 // test_exponential_dist.cpp
  10
  11 #include <boost/math/tools/test.hpp>
  12 #include <boost/math/concepts/real_concept.hpp> // for real_concept
  13 #include <boost/math/distributions/exponential.hpp>
  14     using boost::math::exponential_distribution;
  15
  16 #define BOOST_TEST_MAIN
  17 #include <boost/test/unit_test.hpp> // Boost.Test
  18 #include <boost/test/tools/floating_point_comparison.hpp>
  19 #include "test_out_of_range.hpp"
  20
  21 #include <iostream>
  22    using std::cout;
  23    using std::endl;
  24    using std::setprecision;
  25
  26 template <class RealType>
  27 void test_spot(RealType l, RealType x, RealType p, RealType q, RealType tolerance)
  28 {
  29    BOOST_CHECK_CLOSE(
  30       ::boost::math::cdf(
  31          exponential_distribution<RealType>(l),
  32          x),
  33          p,
  34          tolerance); // %
  35    BOOST_CHECK_CLOSE(
  36       ::boost::math::cdf(
  37          complement(exponential_distribution<RealType>(l),
  38          x)),
  39          q,
  40          tolerance); // %
  41    if(p < 0.999)
  42    {
  43       BOOST_CHECK_CLOSE(
  44          ::boost::math::quantile(
  45             exponential_distribution<RealType>(l),
  46             p),
  47             x,
  48             tolerance); // %
  49    }
  50    if(q < 0.999)
  51    {
  52       BOOST_CHECK_CLOSE(
  53          ::boost::math::quantile(
  54             complement(exponential_distribution<RealType>(l),
  55             q)),
  56             x,
  57             tolerance); // %
  58    }
  59 }
  60
  61 template <class RealType>
  62 void test_spots(RealType T)
  63 {
  64    // Basic sanity checks.
  65    // 50 eps as a percentage, up to a maximum of double precision
  66    // (that's the limit of our test data: obtained by punching
  67    // numbers into a calculator).
  68    RealType tolerance = (std::max)(
  69       static_cast<RealType>(boost::math::tools::epsilon<double>()),
  70       boost::math::tools::epsilon<RealType>());
  71    tolerance *= 50 * 100;
  72    // #  pragma warning(disable: 4100) // unreferenced formal parameter.
  73    // prevent his spurious warning.
  74    if (T != 0)
  75    {
  76      cout << "Expect parameter T == 0!" << endl;
  77    }
  78     cout << "Tolerance for type " << typeid(T).name()  << " is " << tolerance << " %" << endl;
  79
  80    test_spot(
  81       static_cast<RealType>(0.5), // lambda
  82       static_cast<RealType>(0.125), // x
  83       static_cast<RealType>(0.060586937186524213880289175377695L), // p
  84       static_cast<RealType>(0.93941306281347578611971082462231L), //q
  85       tolerance);
  86    test_spot(
  87       static_cast<RealType>(0.5), // lambda
  88       static_cast<RealType>(5), // x
  89       static_cast<RealType>(0.91791500137610120483047132553284L), // p
  90       static_cast<RealType>(0.08208499862389879516952867446716L), //q
  91       tolerance);
  92    test_spot(
  93       static_cast<RealType>(2), // lambda
  94       static_cast<RealType>(0.125), // x
  95       static_cast<RealType>(0.22119921692859513175482973302168L), // p
  96       static_cast<RealType>(0.77880078307140486824517026697832L), //q
  97       tolerance);
  98    test_spot(
  99       static_cast<RealType>(2), // lambda
 100       static_cast<RealType>(5), // x
 101       static_cast<RealType>(0.99995460007023751514846440848444L), // p
 102       static_cast<RealType>(4.5399929762484851535591515560551e-5L), //q
 103       tolerance);
 104
 105    //
 106    // Some spot tests generated by MathCAD pexp(x,r):
 107    //
 108    test_spot(
 109       static_cast<RealType>(1), // lambda
 110       static_cast<RealType>(1), // x
 111       static_cast<RealType>(6.321205588285580E-001L), // p
 112       static_cast<RealType>(1-6.321205588285580E-001L), //q
 113       tolerance);
 114    test_spot(
 115       static_cast<RealType>(2), // lambda
 116       static_cast<RealType>(1), // x
 117       static_cast<RealType>(8.646647167633870E-001L), // p
 118       static_cast<RealType>(1-8.646647167633870E-001L), //q
 119       tolerance);
 120    test_spot(
 121       static_cast<RealType>(1), // lambda
 122       static_cast<RealType>(0.5), // x
 123       static_cast<RealType>(3.934693402873670E-001L), // p
 124       static_cast<RealType>(1-3.934693402873670E-001L), //q
 125       tolerance);
 126    test_spot(
 127       static_cast<RealType>(0.1), // lambda
 128       static_cast<RealType>(1), // x
 129       static_cast<RealType>(9.516258196404040E-002L), // p
 130       static_cast<RealType>(1-9.516258196404040E-002L), //q
 131       tolerance);
 132    test_spot(
 133       static_cast<RealType>(10), // lambda
 134       static_cast<RealType>(1), // x
 135       static_cast<RealType>(9.999546000702380E-001L), // p
 136       static_cast<RealType>(1-9.999546000702380E-001L), //q
 137       tolerance*10000); // we loose four digits to cancellation
 138    test_spot(
 139       static_cast<RealType>(0.1), // lambda
 140       static_cast<RealType>(10), // x
 141       static_cast<RealType>(6.321205588285580E-001L), // p
 142       static_cast<RealType>(1-6.321205588285580E-001L), //q
 143       tolerance);
 144    test_spot(
 145       static_cast<RealType>(1), // lambda
 146       static_cast<RealType>(0.01), // x
 147       static_cast<RealType>(9.950166250831950E-003L), // p
 148       static_cast<RealType>(1-9.950166250831950E-003L), //q
 149       tolerance);
 150    test_spot(
 151       static_cast<RealType>(1), // lambda
 152       static_cast<RealType>(0.0001), // x
 153       static_cast<RealType>(9.999500016666250E-005L), // p
 154       static_cast<RealType>(1-9.999500016666250E-005L), //q
 155       tolerance);
 156    /*
 157    // This test data appears to be erroneous, MathCad appears
 158    // to suffer from cancellation error as x -> 0
 159    test_spot(
 160       static_cast<RealType>(1), // lambda
 161       static_cast<RealType>(0.0000001), // x
 162       static_cast<RealType>(9.999999499998730E-008L), // p
 163       static_cast<RealType>(1-9.999999499998730E-008L), //q
 164       tolerance);
 165    */
 166
 167    BOOST_CHECK_CLOSE(
 168       ::boost::math::pdf(
 169          exponential_distribution<RealType>(0.5),
 170          static_cast<RealType>(0.125)),              // x
 171          static_cast<RealType>(0.46970653140673789305985541231115L),                // probability.
 172          tolerance); // %
 173    BOOST_CHECK_CLOSE(
 174       ::boost::math::pdf(
 175          exponential_distribution<RealType>(0.5),
 176          static_cast<RealType>(5)),              // x
 177          static_cast<RealType>(0.04104249931194939758476433723358L),                // probability.
 178          tolerance); // %
 179    BOOST_CHECK_CLOSE(
 180       ::boost::math::pdf(
 181          exponential_distribution<RealType>(2),
 182          static_cast<RealType>(0.125)),              // x
 183          static_cast<RealType>(1.5576015661428097364903405339566L),                // probability.
 184          tolerance); // %
 185    BOOST_CHECK_CLOSE(
 186       ::boost::math::pdf(
 187          exponential_distribution<RealType>(2),
 188          static_cast<RealType>(5)),              // x
 189          static_cast<RealType>(9.0799859524969703071183031121101e-5L),                // probability.
 190          tolerance); // %
 191
 192    BOOST_CHECK_CLOSE(
 193       ::boost::math::mean(
 194          exponential_distribution<RealType>(2)),
 195          static_cast<RealType>(0.5),
 196          tolerance); // %
 197    BOOST_CHECK_CLOSE(
 198       ::boost::math::standard_deviation(
 199          exponential_distribution<RealType>(2)),
 200          static_cast<RealType>(0.5),
 201          tolerance); // %
 202    BOOST_CHECK_CLOSE(
 203       ::boost::math::mode(
 204          exponential_distribution<RealType>(2)),
 205          static_cast<RealType>(0),
 206          tolerance); // %
 207
 208    BOOST_CHECK_CLOSE(
 209       ::boost::math::median(
 210          exponential_distribution<RealType>(4)),
 211          static_cast<RealType>(0.693147180559945309417232121458176568075500134360255254) / 4,
 212          tolerance); // %
 213
 214    BOOST_CHECK_CLOSE(
 215       ::boost::math::skewness(
 216          exponential_distribution<RealType>(2)),
 217          static_cast<RealType>(2),
 218          tolerance); // %
 219    BOOST_CHECK_CLOSE(
 220       ::boost::math::kurtosis(
 221          exponential_distribution<RealType>(2)),
 222          static_cast<RealType>(9),
 223          tolerance); // %
 224    BOOST_CHECK_CLOSE(
 225       ::boost::math::kurtosis_excess(
 226          exponential_distribution<RealType>(2)),
 227          static_cast<RealType>(6),
 228          tolerance); // %
 229
 230    //
 231    // Things that are errors:
 232    //
 233    exponential_distribution<RealType> dist(0.5);
 234    BOOST_MATH_CHECK_THROW(
 235        quantile(dist, RealType(1.0)),
 236        std::overflow_error);
 237    BOOST_MATH_CHECK_THROW(
 238        quantile(complement(dist, RealType(0.0))),
 239        std::overflow_error);
 240    BOOST_MATH_CHECK_THROW(
 241        pdf(dist, RealType(-1)),
 242        std::domain_error);
 243    BOOST_MATH_CHECK_THROW(
 244        cdf(dist, RealType(-1)),
 245        std::domain_error);
 246    BOOST_MATH_CHECK_THROW(
 247        cdf(exponential_distribution<RealType>(-1), RealType(1)),
 248        std::domain_error);
 249    BOOST_MATH_CHECK_THROW(
 250        quantile(dist, RealType(-1)),
 251        std::domain_error);
 252    BOOST_MATH_CHECK_THROW(
 253        quantile(dist, RealType(2)),
 254        std::domain_error);
 255
 256    check_out_of_range<exponential_distribution<RealType> >(2);
 257    BOOST_MATH_CHECK_THROW(exponential_distribution<RealType>(0), std::domain_error);
 258    BOOST_MATH_CHECK_THROW(exponential_distribution<RealType>(-1), std::domain_error);
 259    if(std::numeric_limits<RealType>::has_infinity)
 260    {
 261       RealType inf = std::numeric_limits<RealType>::infinity();
 262       BOOST_CHECK_EQUAL(pdf(exponential_distribution<RealType>(2), inf), 0);
 263       BOOST_CHECK_EQUAL(cdf(exponential_distribution<RealType>(2), inf), 1);
 264       BOOST_CHECK_EQUAL(cdf(complement(exponential_distribution<RealType>(2), inf)), 0);
 265    }
 266 } // template <class RealType>void test_spots(RealType)
 267
 268 BOOST_AUTO_TEST_CASE( test_main )
 269 {
 270   // Check that can generate exponential distribution using the two convenience methods:
 271    boost::math::exponential mycexp1(1.); // Using typedef
 272    exponential_distribution<> myexp2(1.); // Using default RealType double.
 273
 274     // Basic sanity-check spot values.
 275    // (Parameter value, arbitrarily zero, only communicates the floating point type).
 276   test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 %
 277   test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 %
 278 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 279   test_spots(0.0L); // Test long double.
 280 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 281   test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
 282 #endif
 283 #else
 284    std::cout << "<note>The long double tests have been disabled on this platform "
 285       "either because the long double overloads of the usual math functions are "
 286       "not available at all, or because they are too inaccurate for these tests "
 287       "to pass.</note>" << std::endl;
 288 #endif
 289
 290 } // BOOST_AUTO_TEST_CASE( test_main )
 291
 292 /*
 293
 294 Output is:
 295
 296 Running 1 test case...
 297 Tolerance for type float is 0.000596046 %
 298 Tolerance for type double is 1.11022e-012 %
 299 Tolerance for type long double is 1.11022e-012 %
 300 Tolerance for type class boost::math::concepts::real_concept is 1.11022e-012 %
 301 *** No errors detected
 302
 303 */