libs/math/test/test_weibull.cpp

   1 // Copyright John Maddock 2006, 2012.
   2 // Copyright Paul A. Bristow 2007, 2012.
   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_weibull.cpp
  10
  11 #ifdef _MSC_VER
  12 #  pragma warning (disable : 4127) //  conditional expression is constant.
  13 #endif
  14
  15
  16 #include <boost/math/concepts/real_concept.hpp> // for real_concept
  17 #define BOOST_TEST_MAIN
  18 #include <boost/test/unit_test.hpp> // Boost.Test
  19 #include <boost/test/tools/floating_point_comparison.hpp>
  20
  21 #include <boost/math/distributions/weibull.hpp>
  22     using boost::math::weibull_distribution;
  23 #include <boost/math/tools/test.hpp>
  24 #include "test_out_of_range.hpp"
  25
  26 #include <iostream>
  27    using std::cout;
  28    using std::endl;
  29    using std::setprecision;
  30 #include <limits>
  31   using std::numeric_limits;
  32
  33 template <class RealType>
  34 void check_weibull(RealType shape, RealType scale, RealType x, RealType p, RealType q, RealType tol)
  35 {
  36    BOOST_CHECK_CLOSE(
  37       ::boost::math::cdf(
  38          weibull_distribution<RealType>(shape, scale),       // distribution.
  39          x),                                            // random variable.
  40          p,                                             // probability.
  41          tol);                                          // %tolerance.
  42    BOOST_CHECK_CLOSE(
  43       ::boost::math::cdf(
  44          complement(
  45             weibull_distribution<RealType>(shape, scale),    // distribution.
  46             x)),                                        // random variable.
  47          q,                                             // probability complement.
  48          tol);                                          // %tolerance.
  49    BOOST_CHECK_CLOSE(
  50       ::boost::math::quantile(
  51          weibull_distribution<RealType>(shape, scale),       // distribution.
  52          p),                                            // probability.
  53          x,                                             // random variable.
  54          tol);                                          // %tolerance.
  55    BOOST_CHECK_CLOSE(
  56       ::boost::math::quantile(
  57          complement(
  58             weibull_distribution<RealType>(shape, scale),    // distribution.
  59             q)),                                        // probability complement.
  60          x,                                             // random variable.
  61          tol);                                          // %tolerance.
  62 }
  63
  64 template <class RealType>
  65 void test_spots(RealType)
  66 {
  67    // Basic sanity checks
  68    //
  69    // These test values were generated for the normal distribution
  70    // using the online calculator at
  71    // http://espse.ed.psu.edu/edpsych/faculty/rhale/hale/507Mat/statlets/free/pdist.htm
  72    //
  73    // Tolerance is just over 5 decimal digits expressed as a persentage:
  74    // that's the limit of the test data.
  75    RealType tolerance = 2e-5f * 100;
  76    cout << "Tolerance for type " << typeid(RealType).name()  << " is " << tolerance << " %" << endl;
  77
  78    using std::exp;
  79
  80    check_weibull(
  81       static_cast<RealType>(0.25),     // shape
  82       static_cast<RealType>(0.5),     // scale
  83       static_cast<RealType>(0.1),     // x
  84       static_cast<RealType>(0.487646),   // p
  85       static_cast<RealType>(1-0.487646),   // q
  86       tolerance);
  87    check_weibull(
  88       static_cast<RealType>(0.25),     // shape
  89       static_cast<RealType>(0.5),     // scale
  90       static_cast<RealType>(0.5),     // x
  91       static_cast<RealType>(1-0.367879),   // p
  92       static_cast<RealType>(0.367879),   // q
  93       tolerance);
  94    check_weibull(
  95       static_cast<RealType>(0.25),     // shape
  96       static_cast<RealType>(0.5),     // scale
  97       static_cast<RealType>(1),     // x
  98       static_cast<RealType>(1-0.304463),   // p
  99       static_cast<RealType>(0.304463),   // q
 100       tolerance);
 101    check_weibull(
 102       static_cast<RealType>(0.25),     // shape
 103       static_cast<RealType>(0.5),     // scale
 104       static_cast<RealType>(2),     // x
 105       static_cast<RealType>(1-0.243117),   // p
 106       static_cast<RealType>(0.243117),   // q
 107       tolerance);
 108    check_weibull(
 109       static_cast<RealType>(0.25),     // shape
 110       static_cast<RealType>(0.5),     // scale
 111       static_cast<RealType>(5),     // x
 112       static_cast<RealType>(1-0.168929),   // p
 113       static_cast<RealType>(0.168929),   // q
 114       tolerance);
 115
 116    check_weibull(
 117       static_cast<RealType>(0.5),     // shape
 118       static_cast<RealType>(2),     // scale
 119       static_cast<RealType>(0.1),     // x
 120       static_cast<RealType>(0.200371),   // p
 121       static_cast<RealType>(1-0.200371),   // q
 122       tolerance);
 123    check_weibull(
 124       static_cast<RealType>(0.5),     // shape
 125       static_cast<RealType>(2),     // scale
 126       static_cast<RealType>(0.5),     // x
 127       static_cast<RealType>(0.393469),   // p
 128       static_cast<RealType>(1-0.393469),   // q
 129       tolerance);
 130    check_weibull(
 131       static_cast<RealType>(0.5),     // shape
 132       static_cast<RealType>(2),     // scale
 133       static_cast<RealType>(1),     // x
 134       static_cast<RealType>(1-0.493069),   // p
 135       static_cast<RealType>(0.493069),   // q
 136       tolerance);
 137    check_weibull(
 138       static_cast<RealType>(0.5),     // shape
 139       static_cast<RealType>(2),     // scale
 140       static_cast<RealType>(2),     // x
 141       static_cast<RealType>(1-0.367879),   // p
 142       static_cast<RealType>(0.367879),   // q
 143       tolerance);
 144    check_weibull(
 145       static_cast<RealType>(0.5),     // shape
 146       static_cast<RealType>(2),     // scale
 147       static_cast<RealType>(5),     // x
 148       static_cast<RealType>(1-0.205741),   // p
 149       static_cast<RealType>(0.205741),   // q
 150       tolerance);
 151
 152    check_weibull(
 153       static_cast<RealType>(2),     // shape
 154       static_cast<RealType>(0.25),     // scale
 155       static_cast<RealType>(0.1),     // x
 156       static_cast<RealType>(0.147856),   // p
 157       static_cast<RealType>(1-0.147856),   // q
 158       tolerance);
 159    check_weibull(
 160       static_cast<RealType>(2),     // shape
 161       static_cast<RealType>(0.25),     // scale
 162       static_cast<RealType>(0.5),     // x
 163       static_cast<RealType>(1-0.018316),   // p
 164       static_cast<RealType>(0.018316),   // q
 165       tolerance);
 166
 167    /*
 168    This test value came from
 169    http://espse.ed.psu.edu/edpsych/faculty/rhale/hale/507Mat/statlets/free/pdist.htm
 170    but appears to be grossly incorrect: certainly it does not agree with the values
 171    I get from pushing numbers into a calculator (0.0001249921878255106610615995196123).
 172    Strangely other test values generated for the same shape and scale parameters do look OK.
 173    check_weibull(
 174       static_cast<RealType>(3),     // shape
 175       static_cast<RealType>(2),     // scale
 176       static_cast<RealType>(0.1),     // x
 177       static_cast<RealType>(1.25E-40),   // p
 178       static_cast<RealType>(1-1.25E-40),   // q
 179       tolerance);
 180       */
 181    check_weibull(
 182       static_cast<RealType>(3),     // shape
 183       static_cast<RealType>(2),     // scale
 184       static_cast<RealType>(0.5),     // x
 185       static_cast<RealType>(0.015504),   // p
 186       static_cast<RealType>(1-0.015504),   // q
 187       tolerance * 10); // few digits in test value
 188    check_weibull(
 189       static_cast<RealType>(3),     // shape
 190       static_cast<RealType>(2),     // scale
 191       static_cast<RealType>(1),     // x
 192       static_cast<RealType>(0.117503),   // p
 193       static_cast<RealType>(1-0.117503),   // q
 194       tolerance);
 195    check_weibull(
 196       static_cast<RealType>(3),     // shape
 197       static_cast<RealType>(2),     // scale
 198       static_cast<RealType>(2),     // x
 199       static_cast<RealType>(1-0.367879),   // p
 200       static_cast<RealType>(0.367879),   // q
 201       tolerance);
 202
 203    //
 204    // Tests for PDF
 205    //
 206    BOOST_CHECK_CLOSE(
 207       pdf(weibull_distribution<RealType>(0.25, 0.5), static_cast<RealType>(0.1)),
 208       static_cast<RealType>(0.856579),
 209       tolerance);
 210    BOOST_CHECK_CLOSE(
 211       pdf(weibull_distribution<RealType>(0.25, 0.5), static_cast<RealType>(0.5)),
 212       static_cast<RealType>(0.183940),
 213       tolerance);
 214    BOOST_CHECK_CLOSE(
 215       pdf(weibull_distribution<RealType>(0.25, 0.5), static_cast<RealType>(5)),
 216       static_cast<RealType>(0.015020),
 217       tolerance * 10); // fewer digits in test value
 218    BOOST_CHECK_CLOSE(
 219       pdf(weibull_distribution<RealType>(0.5, 2), static_cast<RealType>(0.1)),
 220       static_cast<RealType>(0.894013),
 221       tolerance);
 222    BOOST_CHECK_CLOSE(
 223       pdf(weibull_distribution<RealType>(0.5, 2), static_cast<RealType>(0.5)),
 224       static_cast<RealType>(0.303265),
 225       tolerance);
 226    BOOST_CHECK_CLOSE(
 227       pdf(weibull_distribution<RealType>(0.5, 2), static_cast<RealType>(1)),
 228       static_cast<RealType>(0.174326),
 229       tolerance);
 230    BOOST_CHECK_CLOSE(
 231       pdf(weibull_distribution<RealType>(2, 0.25), static_cast<RealType>(0.1)),
 232       static_cast<RealType>(2.726860),
 233       tolerance);
 234    BOOST_CHECK_CLOSE(
 235       pdf(weibull_distribution<RealType>(2, 0.25), static_cast<RealType>(0.5)),
 236       static_cast<RealType>(0.293050),
 237       tolerance);
 238    BOOST_CHECK_CLOSE(
 239       pdf(weibull_distribution<RealType>(3, 2), static_cast<RealType>(1)),
 240       static_cast<RealType>(0.330936),
 241       tolerance);
 242    BOOST_CHECK_CLOSE(
 243       pdf(weibull_distribution<RealType>(3, 2), static_cast<RealType>(2)),
 244       static_cast<RealType>(0.551819),
 245       tolerance);
 246
 247    //
 248    // These test values were obtained using the formulas at
 249    // http://en.wikipedia.org/wiki/Weibull_distribution
 250    // which are subtly different to (though mathematically
 251    // the same as) the ones on the Mathworld site
 252    // http://mathworld.wolfram.com/WeibullDistribution.html
 253    // which are the ones used in the implementation.
 254    // The assumption is that if both computation methods
 255    // agree then the implementation is probably correct...
 256    // What's not clear is which method is more accurate.
 257    //
 258    tolerance = (std::max)(
 259       boost::math::tools::epsilon<RealType>(),
 260       static_cast<RealType>(boost::math::tools::epsilon<double>())) * 5 * 100; // 5 eps as a percentage
 261    cout << "Tolerance for type " << typeid(RealType).name()  << " is " << tolerance << " %" << endl;
 262    weibull_distribution<RealType> dist(2, 3);
 263    RealType x = static_cast<RealType>(0.125);
 264
 265    BOOST_MATH_STD_USING // ADL of std lib math functions
 266
 267    // mean:
 268    BOOST_CHECK_CLOSE(
 269       mean(dist)
 270       , dist.scale() * boost::math::tgamma(1 + 1 / dist.shape()), tolerance);
 271    // variance:
 272    BOOST_CHECK_CLOSE(
 273       variance(dist)
 274       , dist.scale() * dist.scale() * boost::math::tgamma(1 + 2 / dist.shape()) - mean(dist) * mean(dist), tolerance);
 275    // std deviation:
 276    BOOST_CHECK_CLOSE(
 277     standard_deviation(dist)
 278     , sqrt(variance(dist)), tolerance);
 279    // hazard:
 280    BOOST_CHECK_CLOSE(
 281     hazard(dist, x)
 282     , pdf(dist, x) / cdf(complement(dist, x)), tolerance);
 283    // cumulative hazard:
 284    BOOST_CHECK_CLOSE(
 285     chf(dist, x)
 286     , -log(cdf(complement(dist, x))), tolerance);
 287    // coefficient_of_variation:
 288    BOOST_CHECK_CLOSE(
 289     coefficient_of_variation(dist)
 290     , standard_deviation(dist) / mean(dist), tolerance);
 291    // mode:
 292    BOOST_CHECK_CLOSE(
 293     mode(dist)
 294     , dist.scale() * pow((dist.shape() - 1) / dist.shape(), 1/dist.shape()), tolerance);
 295    // median:
 296    BOOST_CHECK_CLOSE(
 297     median(dist)
 298     , dist.scale() * pow(log(static_cast<RealType>(2)), 1 / dist.shape()), tolerance);
 299    // skewness:
 300    BOOST_CHECK_CLOSE(
 301     skewness(dist),
 302     (boost::math::tgamma(1 + 3/dist.shape()) * pow(dist.scale(), RealType(3)) - 3 * mean(dist) * variance(dist) - pow(mean(dist), RealType(3))) / pow(standard_deviation(dist), RealType(3)),
 303     tolerance * 100);
 304    // kertosis:
 305    BOOST_CHECK_CLOSE(
 306     kurtosis(dist)
 307     , kurtosis_excess(dist) + 3, tolerance);
 308    // kertosis excess:
 309    BOOST_CHECK_CLOSE(
 310     kurtosis_excess(dist),
 311     (pow(dist.scale(), RealType(4)) * boost::math::tgamma(1 + 4/dist.shape())
 312          - 3 * variance(dist) * variance(dist)
 313          - 4 * skewness(dist) * variance(dist) * standard_deviation(dist) * mean(dist)
 314          - 6 * variance(dist) * mean(dist) * mean(dist)
 315          - pow(mean(dist), RealType(4))) / (variance(dist) * variance(dist)),
 316     tolerance * 1000);
 317
 318    //
 319    // Special cases:
 320    //
 321    BOOST_CHECK(cdf(dist, 0) == 0);
 322    BOOST_CHECK(cdf(complement(dist, 0)) == 1);
 323    BOOST_CHECK(quantile(dist, 0) == 0);
 324    BOOST_CHECK(quantile(complement(dist, 1)) == 0);
 325
 326    BOOST_CHECK_EQUAL(pdf(weibull_distribution<RealType>(1, 1), 0), 1);
 327
 328    //
 329    // Error checks:
 330    //
 331    BOOST_MATH_CHECK_THROW(weibull_distribution<RealType>(1, -1), std::domain_error);
 332    BOOST_MATH_CHECK_THROW(weibull_distribution<RealType>(-1, 1), std::domain_error);
 333    BOOST_MATH_CHECK_THROW(weibull_distribution<RealType>(1, 0), std::domain_error);
 334    BOOST_MATH_CHECK_THROW(weibull_distribution<RealType>(0, 1), std::domain_error);
 335    BOOST_MATH_CHECK_THROW(pdf(dist, -1), std::domain_error);
 336    BOOST_MATH_CHECK_THROW(cdf(dist, -1), std::domain_error);
 337    BOOST_MATH_CHECK_THROW(cdf(complement(dist, -1)), std::domain_error);
 338    BOOST_MATH_CHECK_THROW(quantile(dist, 1), std::overflow_error);
 339    BOOST_MATH_CHECK_THROW(quantile(complement(dist, 0)), std::overflow_error);
 340    BOOST_MATH_CHECK_THROW(quantile(dist, -1), std::domain_error);
 341    BOOST_MATH_CHECK_THROW(quantile(complement(dist, -1)), std::domain_error);
 342
 343    BOOST_CHECK_EQUAL(pdf(dist, 0), exp(-pow(RealType(0) / RealType(3), RealType(2))) * pow(RealType(0), RealType(1)) * RealType(2) / RealType(3));
 344    BOOST_CHECK_EQUAL(pdf(weibull_distribution<RealType>(1, 3), 0), exp(-pow(RealType(0) / RealType(3), RealType(1))) * pow(RealType(0), RealType(0)) * RealType(1) / RealType(3));
 345    BOOST_MATH_CHECK_THROW(pdf(weibull_distribution<RealType>(0.5, 3), 0), std::overflow_error);
 346
 347    check_out_of_range<weibull_distribution<RealType> >(1, 1);
 348 } // template <class RealType>void test_spots(RealType)
 349
 350 BOOST_AUTO_TEST_CASE( test_main )
 351 {
 352
 353   // Check that can construct weibull distribution using the two convenience methods:
 354   using namespace boost::math;
 355   weibull myw1(2); // Using typedef
 356    weibull_distribution<> myw2(2); // Using default RealType double.
 357
 358     // Basic sanity-check spot values.
 359    // (Parameter value, arbitrarily zero, only communicates the floating point type).
 360   test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 %
 361   test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 %
 362 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 363   test_spots(0.0L); // Test long double.
 364 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x0582))
 365   test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
 366 #endif
 367 #else
 368    std::cout << "<note>The long double tests have been disabled on this platform "
 369       "either because the long double overloads of the usual math functions are "
 370       "not available at all, or because they are too inaccurate for these tests "
 371       "to pass.</note>" << std::endl;
 372 #endif
 373
 374
 375 } // BOOST_AUTO_TEST_CASE( test_main )
 376
 377 /*
 378
 379 Output:
 380
 381   Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Debug\test_weibull.exe"
 382   Running 1 test case...
 383   Tolerance for type float is 0.002 %
 384   Tolerance for type float is 5.96046e-005 %
 385   Tolerance for type double is 0.002 %
 386   Tolerance for type double is 1.11022e-013 %
 387   Tolerance for type long double is 0.002 %
 388   Tolerance for type long double is 1.11022e-013 %
 389   Tolerance for type class boost::math::concepts::real_concept is 0.002 %
 390   Tolerance for type class boost::math::concepts::real_concept is 1.11022e-013 %
 391
 392   *** No errors detected
 393
 394
 395 */
 396
 397
 398
 399