libs/math/test/test_pareto.cpp

   1 // Copyright Paul A. Bristow 2007, 2009.
   2 // Copyright John Maddock 2006.
   3 // Use, modification and distribution are subject to the
   4 // Boost Software License, Version 1.0.
   5 // (See accompanying file LICENSE_1_0.txt
   6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 // test_pareto.cpp
   9
  10 // http://en.wikipedia.org/wiki/pareto_distribution
  11 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm
  12 // Also:
  13 // Weisstein, Eric W. "pareto Distribution."
  14 // From MathWorld--A Wolfram Web Resource.
  15 // http://mathworld.wolfram.com/paretoDistribution.html
  16
  17
  18 #ifdef _MSC_VER
  19 #  pragma warning(disable: 4127) // conditional expression is constant.
  20 # pragma warning (disable : 4996) // POSIX name for this item is deprecated
  21 # pragma warning (disable : 4224) // nonstandard extension used : formal parameter 'arg' was previously defined as a type
  22 # pragma warning (disable : 4180) // qualifier applied to function type has no meaning; ignored
  23 #  pragma warning(disable: 4100) // unreferenced formal parameter.
  24 #endif
  25
  26 #include <boost/math/tools/test.hpp> // for real_concept
  27 #include <boost/math/concepts/real_concept.hpp> // for real_concept
  28 #define BOOST_TEST_MAIN
  29 #include <boost/test/unit_test.hpp> // Boost.Test
  30 #include <boost/test/tools/floating_point_comparison.hpp>
  31
  32 #include <boost/math/distributions/pareto.hpp>
  33     using boost::math::pareto_distribution;
  34 #include <boost/math/tools/test.hpp>
  35 #include "test_out_of_range.hpp"
  36
  37 #include <iostream>
  38    using std::cout;
  39    using std::endl;
  40    using std::setprecision;
  41 #include <limits>
  42   using std::numeric_limits;
  43
  44   template <class RealType>
  45   void check_pareto(RealType scale, RealType shape, RealType x, RealType p, RealType q, RealType tol)
  46   {
  47     BOOST_CHECK_CLOSE_FRACTION(
  48       ::boost::math::cdf(
  49       pareto_distribution<RealType>(scale, shape),   // distribution.
  50       x),                                            // random variable.
  51       p,                                             // probability.
  52       tol);                                          // tolerance eps.
  53     BOOST_CHECK_CLOSE_FRACTION(
  54       ::boost::math::cdf(
  55       complement(
  56       pareto_distribution<RealType>(scale, shape),   // distribution.
  57       x)),                                           // random variable.
  58       q,                                             // probability complement.
  59       tol);                                          // tolerance eps.
  60     BOOST_CHECK_CLOSE_FRACTION(
  61       ::boost::math::quantile(
  62       pareto_distribution<RealType>(scale, shape),   // distribution.
  63       p),                                            // probability.
  64       x,                                             // random variable.
  65       tol);                                          // tolerance eps.
  66     BOOST_CHECK_CLOSE_FRACTION(
  67       ::boost::math::quantile(
  68       complement(
  69       pareto_distribution<RealType>(scale, shape),    // distribution.
  70       q)),                                        // probability complement.
  71       x,                                             // random variable.
  72       tol);                                          // tolerance eps.
  73   } // check_pareto
  74
  75 template <class RealType>
  76 void test_spots(RealType)
  77 {
  78    // Basic sanity checks.
  79    //
  80    // Tolerance are based on units of epsilon, but capped at
  81    // double precision, since that's the limit of our test data:
  82    //
  83    RealType tol = (std::max)((RealType)boost::math::tools::epsilon<double>(), boost::math::tools::epsilon<RealType>());
  84    RealType tol5eps = tol * 5;
  85    RealType tol10eps = tol * 10;
  86    RealType tol100eps = tol * 100;
  87    RealType tol1000eps = tol * 1000;
  88
  89    check_pareto(
  90       static_cast<RealType>(1.1L), //
  91       static_cast<RealType>(5.5L),
  92       static_cast<RealType>(2.2L),
  93       static_cast<RealType>(0.97790291308792L),
  94       static_cast<RealType>(0.0220970869120796L),
  95       tol10eps * 4);
  96
  97    check_pareto(
  98       static_cast<RealType>(0.5L),
  99       static_cast<RealType>(10.1L),
 100       static_cast<RealType>(1.5L),
 101       static_cast<RealType>(0.99998482686481L),
 102       static_cast<RealType>(1.51731351900608e-005L),
 103       tol100eps * 1000); // Much less accurate as p close to unity.
 104
 105    check_pareto(
 106       static_cast<RealType>(0.1L),
 107       static_cast<RealType>(2.3L),
 108       static_cast<RealType>(1.5L),
 109       static_cast<RealType>(0.99802762220697L),
 110       static_cast<RealType>(0.00197237779302972L),
 111       tol1000eps);
 112
 113    // Example from 23.3 page 259
 114    check_pareto(
 115       static_cast<RealType>(2.30444301457005L),
 116       static_cast<RealType>(4),
 117       static_cast<RealType>(2.4L),
 118       static_cast<RealType>(0.15L),
 119       static_cast<RealType>(0.85L),
 120       tol100eps);
 121
 122    check_pareto(
 123       static_cast<RealType>(2),
 124       static_cast<RealType>(3),
 125       static_cast<RealType>(3.4L),
 126       static_cast<RealType>(0.796458375737838L),
 127       static_cast<RealType>(0.203541624262162L),
 128       tol10eps);
 129
 130    check_pareto( // Probability near 0.5
 131       static_cast<RealType>(2),
 132       static_cast<RealType>(2),
 133       static_cast<RealType>(3),
 134       static_cast<RealType>(0.5555555555555555555555555555555555555556L),
 135       static_cast<RealType>(0.4444444444444444444444444444444444444444L),
 136       tol5eps); // accurate.
 137
 138
 139    // Tests for:
 140
 141    // pdf for shapes 1, 2 & 3 (exact)
 142    BOOST_CHECK_CLOSE_FRACTION(
 143       pdf(pareto_distribution<RealType>(1, 1), 1),
 144       static_cast<RealType>(1), //
 145       tol5eps);
 146
 147     BOOST_CHECK_CLOSE_FRACTION(   pdf(pareto_distribution<RealType>(1, 2), 1),
 148       static_cast<RealType>(2), //
 149       tol5eps);
 150
 151      BOOST_CHECK_CLOSE_FRACTION(   pdf(pareto_distribution<RealType>(1, 3), 1),
 152       static_cast<RealType>(3), //
 153       tol5eps);
 154
 155    // cdf
 156    BOOST_CHECK_EQUAL( // x = scale
 157       cdf(pareto_distribution<RealType>(1, 1), 1),
 158       static_cast<RealType>(0) );
 159
 160    // Compare with values from StatCalc K. Krishnamoorthy,  ISBN 1-58488-635-8 eq 23.1.3
 161    BOOST_CHECK_CLOSE_FRACTION( // small x
 162       cdf(pareto_distribution<RealType>(2, 5), static_cast<RealType>(3.4)),
 163       static_cast<RealType>(0.929570372227626L), tol5eps);
 164
 165    BOOST_CHECK_CLOSE_FRACTION( // small x
 166       cdf(pareto_distribution<RealType>(2, 5), static_cast<RealType>(3.4)),
 167       static_cast<RealType>(1 - 0.0704296277723743L), tol5eps);
 168
 169    BOOST_CHECK_CLOSE_FRACTION( // small x
 170       cdf(complement(pareto_distribution<RealType>(2, 5), static_cast<RealType>(3.4))),
 171       static_cast<RealType>(0.0704296277723743L), tol5eps);
 172
 173    // quantile
 174    BOOST_CHECK_EQUAL( // x = scale
 175       quantile(pareto_distribution<RealType>(1, 1), 0),
 176       static_cast<RealType>(1) );
 177
 178    BOOST_CHECK_EQUAL( // x = scale
 179       quantile(complement(pareto_distribution<RealType>(1, 1), 1)),
 180       static_cast<RealType>(1) );
 181
 182    BOOST_CHECK_CLOSE_FRACTION( // small x
 183       cdf(complement(pareto_distribution<RealType>(2, 5), static_cast<RealType>(3.4))),
 184       static_cast<RealType>(0.0704296277723743L), tol5eps);
 185
 186     using namespace std; // ADL of std names.
 187
 188     pareto_distribution<RealType> pareto15(1, 5);
 189     // Note: shape must be big enough (5) that all moments up to kurtosis are defined
 190     // to allow all functions to be tested.
 191
 192     // mean:
 193     BOOST_CHECK_CLOSE_FRACTION(
 194        mean(pareto15), static_cast<RealType>(1.25), tol5eps); // 1.25 == 5/4
 195     BOOST_CHECK_EQUAL(
 196        mean(pareto15), static_cast<RealType>(1.25)); // 1.25 == 5/4 (expect exact so check equal)
 197
 198     pareto_distribution<RealType> p12(1, 2); //
 199     BOOST_CHECK_EQUAL(
 200        mean(p12), static_cast<RealType>(2)); // Exactly two.
 201
 202     // variance:
 203    BOOST_CHECK_CLOSE_FRACTION(
 204        variance(pareto15), static_cast<RealType>(0.10416666666666667L), tol5eps);
 205     // std deviation:
 206     BOOST_CHECK_CLOSE_FRACTION(
 207        standard_deviation(pareto15), static_cast<RealType>(0.32274861218395140L), tol5eps);
 208     // hazard:   No independent test values found yet.
 209     //BOOST_CHECK_CLOSE_FRACTION(
 210     //   hazard(pareto15, x), pdf(pareto15, x) / cdf(complement(pareto15, x)), tol5eps);
 211     //// cumulative hazard:
 212     //BOOST_CHECK_CLOSE_FRACTION(
 213     //   chf(pareto15, x), -log(cdf(complement(pareto15, x))), tol5eps);
 214     //// coefficient_of_variation:
 215     BOOST_CHECK_CLOSE_FRACTION(
 216        coefficient_of_variation(pareto15), static_cast<RealType>(0.25819888974716110L), tol5eps);
 217     // mode:
 218     BOOST_CHECK_CLOSE_FRACTION(
 219        mode(pareto15), static_cast<RealType>(1), tol5eps);
 220
 221     BOOST_CHECK_CLOSE_FRACTION(
 222        median(pareto15), static_cast<RealType>(1.1486983549970351L), tol5eps);
 223
 224     // skewness:
 225     BOOST_CHECK_CLOSE_FRACTION(
 226        skewness(pareto15), static_cast<RealType>(4.6475800154489004L), tol5eps);
 227     // kertosis:
 228     BOOST_CHECK_CLOSE_FRACTION(
 229        kurtosis(pareto15), static_cast<RealType>(73.8L), tol5eps);
 230     // kertosis excess:
 231     BOOST_CHECK_CLOSE_FRACTION(
 232        kurtosis_excess(pareto15), static_cast<RealType>(70.8L), tol5eps);
 233     // Check difference between kurtosis and excess:
 234     BOOST_CHECK_CLOSE_FRACTION(
 235       kurtosis_excess(pareto15), kurtosis(pareto15) - static_cast<RealType>(3L), tol5eps);
 236     // Check kurtosis excess = kurtosis - 3;
 237
 238     // Error condition checks:
 239     check_out_of_range<pareto_distribution<RealType> >(1, 1);
 240     BOOST_MATH_CHECK_THROW(pdf(pareto_distribution<RealType>(0, 1), 0), std::domain_error);
 241     BOOST_MATH_CHECK_THROW(pdf(pareto_distribution<RealType>(1, 0), 0), std::domain_error);
 242     BOOST_MATH_CHECK_THROW(pdf(pareto_distribution<RealType>(-1, 1), 0), std::domain_error);
 243     BOOST_MATH_CHECK_THROW(pdf(pareto_distribution<RealType>(1, -1), 0), std::domain_error);
 244
 245     BOOST_MATH_CHECK_THROW(pdf(pareto_distribution<RealType>(1, 1), 0), std::domain_error);
 246     BOOST_MATH_CHECK_THROW(cdf(pareto_distribution<RealType>(1, 1), 0), std::domain_error);
 247
 248     BOOST_MATH_CHECK_THROW(quantile(pareto_distribution<RealType>(1, 1), -1), std::domain_error);
 249     BOOST_MATH_CHECK_THROW(quantile(pareto_distribution<RealType>(1, 1), 2), std::domain_error);
 250 } // template <class RealType>void test_spots(RealType)
 251
 252 BOOST_AUTO_TEST_CASE( test_main )
 253 {
 254   // Check that can generate pareto distribution using the two convenience methods:
 255    boost::math::pareto myp1(1., 1); // Using typedef
 256    pareto_distribution<> myp2(1., 1); // Using default RealType double.
 257   boost::math::pareto pareto11; // Use default values (scale = 1, shape = 1).
 258   // Note NOT pareto11() as the compiler will interpret as a function!
 259    // Basic sanity-check spot values.
 260
 261   BOOST_CHECK_EQUAL(pareto11.scale(), 1); // Check defaults again.
 262   BOOST_CHECK_EQUAL(pareto11.shape(), 1);
 263   BOOST_CHECK_EQUAL(myp1.scale(), 1);
 264   BOOST_CHECK_EQUAL(myp1.shape(), 1);
 265   BOOST_CHECK_EQUAL(myp2.scale(), 1);
 266   BOOST_CHECK_EQUAL(myp2.shape(), 1);
 267
 268   // Test range and support using double only,
 269   // because it supports numeric_limits max for pseudo-infinity.
 270   BOOST_CHECK_EQUAL(range(myp2).first, 0); // range 0 to +infinity
 271   BOOST_CHECK_EQUAL(range(myp2).second, (numeric_limits<double>::max)());
 272   BOOST_CHECK_EQUAL(support(myp2).first, myp2.scale()); // support scale to + infinity.
 273   BOOST_CHECK_EQUAL(support(myp2).second, (numeric_limits<double>::max)());
 274
 275   // Check some bad parameters to the distribution.
 276 #ifndef BOOST_NO_EXCEPTIONS
 277    BOOST_MATH_CHECK_THROW(boost::math::pareto mypm1(-1, 1), std::domain_error); // Using typedef
 278    BOOST_MATH_CHECK_THROW(boost::math::pareto myp0(0, 1), std::domain_error); // Using typedef
 279    BOOST_MATH_CHECK_THROW(boost::math::pareto myp1m1(1, -1), std::domain_error); // Using typedef
 280    BOOST_MATH_CHECK_THROW(boost::math::pareto myp10(1, 0), std::domain_error); // Using typedef
 281 #else
 282    BOOST_MATH_CHECK_THROW(boost::math::pareto(-1, 1), std::domain_error); // Using typedef
 283    BOOST_MATH_CHECK_THROW(boost::math::pareto(0, 1), std::domain_error); // Using typedef
 284    BOOST_MATH_CHECK_THROW(boost::math::pareto(1, -1), std::domain_error); // Using typedef
 285    BOOST_MATH_CHECK_THROW(boost::math::pareto(1, 0), std::domain_error); // Using typedef
 286 #endif
 287
 288   // Check some moments that should fail because shape not big enough.
 289   BOOST_MATH_CHECK_THROW(variance(myp2), std::domain_error);
 290   BOOST_MATH_CHECK_THROW(standard_deviation(myp2), std::domain_error);
 291   BOOST_MATH_CHECK_THROW(skewness(myp2), std::domain_error);
 292   BOOST_MATH_CHECK_THROW(kurtosis(myp2), std::domain_error);
 293   BOOST_MATH_CHECK_THROW(kurtosis_excess(myp2), std::domain_error);
 294
 295   // Test on extreme values of distribution parameters,
 296   // using just double because it has numeric_limit infinity etc.
 297 #ifndef BOOST_NO_EXCEPTIONS
 298   BOOST_MATH_CHECK_THROW(boost::math::pareto mypinf1(+std::numeric_limits<double>::infinity(), 1), std::domain_error); // Using typedef
 299    BOOST_MATH_CHECK_THROW(boost::math::pareto myp1inf(1, +std::numeric_limits<double>::infinity()), std::domain_error); // Using typedef
 300    BOOST_MATH_CHECK_THROW(boost::math::pareto mypinf1(+std::numeric_limits<double>::infinity(), +std::numeric_limits<double>::infinity()), std::domain_error); // Using typedef
 301 #else
 302   BOOST_MATH_CHECK_THROW(boost::math::pareto(+std::numeric_limits<double>::infinity(), 1), std::domain_error); // Using typedef
 303    BOOST_MATH_CHECK_THROW(boost::math::pareto(1, +std::numeric_limits<double>::infinity()), std::domain_error); // Using typedef
 304    BOOST_MATH_CHECK_THROW(boost::math::pareto(+std::numeric_limits<double>::infinity(), +std::numeric_limits<double>::infinity()), std::domain_error); // Using typedef
 305 #endif
 306
 307   // Test on extreme values of random variate x, using just double because it has numeric_limit infinity etc..
 308   // No longer allow x to be + or - infinity, then these tests should throw.
 309   BOOST_MATH_CHECK_THROW(pdf(pareto11, +std::numeric_limits<double>::infinity()), std::domain_error); // x = + infinity
 310   BOOST_MATH_CHECK_THROW(pdf(pareto11, -std::numeric_limits<double>::infinity()), std::domain_error); // x = - infinity
 311   BOOST_MATH_CHECK_THROW(cdf(pareto11, +std::numeric_limits<double>::infinity()), std::domain_error); // x = + infinity
 312   BOOST_MATH_CHECK_THROW(cdf(pareto11, -std::numeric_limits<double>::infinity()), std::domain_error); // x = - infinity
 313
 314   BOOST_CHECK_EQUAL(pdf(pareto11, 0.5), 0); // x < scale but > 0
 315   BOOST_CHECK_EQUAL(pdf(pareto11, (std::numeric_limits<double>::min)()), 0); // x almost zero but > 0
 316   BOOST_CHECK_EQUAL(pdf(pareto11, 1), 1); // x == scale, result == shape == 1
 317   BOOST_CHECK_EQUAL(pdf(pareto11, +(std::numeric_limits<double>::max)()), 0); // x = +max, pdf has fallen to zero.
 318
 319   BOOST_MATH_CHECK_THROW(pdf(pareto11, 0), std::domain_error); // x == 0
 320   BOOST_MATH_CHECK_THROW(pdf(pareto11, -1), std::domain_error); // x = -1
 321   BOOST_MATH_CHECK_THROW(pdf(pareto11, -(std::numeric_limits<double>::max)()), std::domain_error); // x = - max
 322   BOOST_MATH_CHECK_THROW(pdf(pareto11, -(std::numeric_limits<double>::min)()), std::domain_error); // x = - min
 323
 324   BOOST_CHECK_EQUAL(cdf(pareto11, 1), 0); // x == scale, cdf = zero.
 325   BOOST_CHECK_EQUAL(cdf(pareto11, +(std::numeric_limits<double>::max)()), 1); // x = + max, cdf = unity.
 326
 327   BOOST_MATH_CHECK_THROW(cdf(pareto11, 0), std::domain_error); // x == 0
 328   BOOST_MATH_CHECK_THROW(cdf(pareto11, -(std::numeric_limits<double>::min)()), std::domain_error); // x = - min,
 329   BOOST_MATH_CHECK_THROW(cdf(pareto11, -(std::numeric_limits<double>::max)()), std::domain_error); // x = - max,
 330
 331    // (Parameter value, arbitrarily zero, only communicates the floating point type).
 332   test_spots(0.0F); // Test float. OK at decdigits = 0 tol5eps = 0.0001 %
 333   test_spots(0.0); // Test double. OK at decdigits 7, tol5eps = 1e07 %
 334 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 335   test_spots(0.0L); // Test long double.
 336 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x0582))
 337   test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
 338 #endif
 339 #else
 340    std::cout << "<note>The long double tests have been disabled on this platform "
 341       "either because the long double overloads of the usual math functions are "
 342       "not available at all, or because they are too inaccurate for these tests "
 343       "to pass.</note>" << std::endl;
 344 #endif
 345
 346
 347 } // BOOST_AUTO_TEST_CASE( test_main )
 348
 349 /*
 350
 351 Output:
 352
 353 Compiling...
 354 test_pareto.cpp
 355 Linking...
 356 Embedding manifest...
 357 Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_pareto.exe"
 358 Running 1 test case...
 359 *** No errors detected
 360
 361
 362
 363 */
 364
 365