ceph/src/boost/libs/math/test/test_lognormal.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_lognormal.cpp
  10
  11 #include <boost/math/concepts/real_concept.hpp> // for real_concept
  12 #define BOOST_TEST_MAIN
  13 #include <boost/test/unit_test.hpp> // Boost.Test
  14 #include <boost/test/floating_point_comparison.hpp>
  15
  16 #include <boost/math/distributions/lognormal.hpp>
  17     using boost::math::lognormal_distribution;
  18 #include <boost/math/tools/test.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 #include <limits>
  26   using std::numeric_limits;
  27 #include <cassert>
  28
  29 template <class RealType>
  30 void check_lognormal(RealType loc, RealType scale, RealType x, RealType p, RealType q, RealType tol)
  31 {
  32    BOOST_CHECK_CLOSE(
  33       ::boost::math::cdf(
  34          lognormal_distribution<RealType>(loc, scale),       // distribution.
  35          x),                                            // random variable.
  36          p,                                             // probability.
  37          tol);                                          // %tolerance.
  38    BOOST_CHECK_CLOSE(
  39       ::boost::math::cdf(
  40          complement(
  41             lognormal_distribution<RealType>(loc, scale),    // distribution.
  42             x)),                                        // random variable.
  43          q,                                             // probability complement.
  44          tol);                                          // %tolerance.
  45    BOOST_CHECK_CLOSE(
  46       ::boost::math::quantile(
  47          lognormal_distribution<RealType>(loc, scale),       // distribution.
  48          p),                                            // probability.
  49          x,                                             // random variable.
  50          tol);                                          // %tolerance.
  51    BOOST_CHECK_CLOSE(
  52       ::boost::math::quantile(
  53          complement(
  54             lognormal_distribution<RealType>(loc, scale),    // distribution.
  55             q)),                                        // probability complement.
  56          x,                                             // random variable.
  57          tol);                                          // %tolerance.
  58 }
  59
  60 template <class RealType>
  61 void test_spots(RealType)
  62 {
  63
  64    // Basic sanity checks.
  65    RealType tolerance = 5e-3 * 100;
  66    // Some tests only pass at 1e-4 because values generated by
  67    // http://faculty.vassar.edu/lowry/VassarStats.html
  68    // give only 5 or 6 *fixed* places, so small values have fewer digits.
  69
  70    cout << "Tolerance for type " << typeid(RealType).name()  << " is " << tolerance << " %" << endl;
  71
  72    using std::exp;
  73
  74    //
  75    // These test values were generated for the normal distribution
  76    // using the online calculator at http://faculty.vassar.edu/lowry/VassarStats.html
  77    // and then exponentiating the random variate x.
  78    //
  79
  80    check_lognormal(
  81       static_cast<RealType>(0),     // location
  82       static_cast<RealType>(5),     // scale
  83       static_cast<RealType>(1),     // x
  84       static_cast<RealType>(0.5),   // p
  85       static_cast<RealType>(0.5),   // q
  86       tolerance);
  87    check_lognormal(
  88       static_cast<RealType>(2),     // location
  89       static_cast<RealType>(2),     // scale
  90       static_cast<RealType>(exp(1.8)),     // x
  91       static_cast<RealType>(0.46017),   // p
  92       static_cast<RealType>(1-0.46017),   // q
  93       tolerance);
  94    check_lognormal(
  95       static_cast<RealType>(2),     // location
  96       static_cast<RealType>(2),     // scale
  97       static_cast<RealType>(exp(2.2)),     // x
  98       static_cast<RealType>(1-0.46017),   // p
  99       static_cast<RealType>(0.46017),   // q
 100       tolerance);
 101    check_lognormal(
 102       static_cast<RealType>(2),     // location
 103       static_cast<RealType>(2),     // scale
 104       static_cast<RealType>(exp(-1.4)),     // x
 105       static_cast<RealType>(0.04457),   // p
 106       static_cast<RealType>(1-0.04457),   // q
 107       tolerance);
 108    check_lognormal(
 109       static_cast<RealType>(2),     // location
 110       static_cast<RealType>(2),     // scale
 111       static_cast<RealType>(exp(5.4)),     // x
 112       static_cast<RealType>(1-0.04457),   // p
 113       static_cast<RealType>(0.04457),   // q
 114       tolerance);
 115
 116    check_lognormal(
 117       static_cast<RealType>(-3),     // location
 118       static_cast<RealType>(5),     // scale
 119       static_cast<RealType>(exp(-5.0)),     // x
 120       static_cast<RealType>(0.34458),   // p
 121       static_cast<RealType>(1-0.34458),   // q
 122       tolerance);
 123    check_lognormal(
 124       static_cast<RealType>(-3),     // location
 125       static_cast<RealType>(5),     // scale
 126       static_cast<RealType>(exp(-1.0)),     // x
 127       static_cast<RealType>(1-0.34458),   // p
 128       static_cast<RealType>(0.34458),   // q
 129       tolerance);
 130    check_lognormal(
 131       static_cast<RealType>(-3),     // location
 132       static_cast<RealType>(5),     // scale
 133       static_cast<RealType>(exp(-9.0)),     // x
 134       static_cast<RealType>(0.11507),   // p
 135       static_cast<RealType>(1-0.11507),   // q
 136       tolerance);
 137    check_lognormal(
 138       static_cast<RealType>(-3),     // location
 139       static_cast<RealType>(5),     // scale
 140       static_cast<RealType>(exp(3.0)),     // x
 141       static_cast<RealType>(1-0.11507),   // p
 142       static_cast<RealType>(0.11507),   // q
 143       tolerance);
 144
 145    //
 146    // Tests for PDF
 147    //
 148    tolerance = boost::math::tools::epsilon<RealType>() * 5 * 100; // 5 eps as a percentage
 149    cout << "Tolerance for type " << typeid(RealType).name()  << " is " << tolerance << " %" << endl;
 150    BOOST_CHECK_CLOSE(
 151       pdf(lognormal_distribution<RealType>(), static_cast<RealType>(1)),
 152       static_cast<RealType>(0.3989422804014326779399460599343818684759L), // 1/sqrt(2*pi)
 153       tolerance);
 154    BOOST_CHECK_CLOSE(
 155       pdf(lognormal_distribution<RealType>(3), exp(static_cast<RealType>(3))),
 156       static_cast<RealType>(0.3989422804014326779399460599343818684759L) / exp(static_cast<RealType>(3)),
 157       tolerance);
 158    BOOST_CHECK_CLOSE(
 159       pdf(lognormal_distribution<RealType>(3, 5), exp(static_cast<RealType>(3))),
 160       static_cast<RealType>(0.3989422804014326779399460599343818684759L / (5 * exp(static_cast<RealType>(3)))),
 161       tolerance);
 162    //
 163    // Spot checks for location = -5, scale = 6,
 164    // use relation to normal to test:
 165    //
 166    for(RealType x = -15; x < 5; x += 0.125)
 167    {
 168       BOOST_CHECK_CLOSE(
 169          pdf(lognormal_distribution<RealType>(-5, 6), exp(x)),
 170          pdf(boost::math::normal_distribution<RealType>(-5, 6), x) / exp(x),
 171          tolerance);
 172    }
 173
 174    //
 175    // These test values were obtained by punching numbers into
 176    // a calculator, using the formulas at http://mathworld.wolfram.com/LogNormalDistribution.html
 177    //
 178    tolerance = (std::max)(
 179       boost::math::tools::epsilon<RealType>(),
 180       static_cast<RealType>(boost::math::tools::epsilon<double>())) * 5 * 100; // 5 eps as a percentage
 181    cout << "Tolerance for type " << typeid(RealType).name()  << " is " << tolerance << " %" << endl;
 182    lognormal_distribution<RealType> dist(8, 3);
 183    RealType x = static_cast<RealType>(0.125);
 184
 185    BOOST_MATH_STD_USING // ADL of std math lib names.
 186
 187    // mean:
 188    BOOST_CHECK_CLOSE(
 189       mean(dist)
 190       , static_cast<RealType>(268337.28652087445695647967378715L), tolerance);
 191    // variance:
 192    BOOST_CHECK_CLOSE(
 193       variance(dist)
 194       , static_cast<RealType>(583389737628117.49553037857325892L), tolerance);
 195    // std deviation:
 196    BOOST_CHECK_CLOSE(
 197     standard_deviation(dist)
 198     , static_cast<RealType>(24153462.228594009489719473727471L), tolerance);
 199    // hazard:
 200    BOOST_CHECK_CLOSE(
 201     hazard(dist, x)
 202     , pdf(dist, x) / cdf(complement(dist, x)), tolerance);
 203    // cumulative hazard:
 204    BOOST_CHECK_CLOSE(
 205     chf(dist, x)
 206     , -log(cdf(complement(dist, x))), tolerance);
 207    // coefficient_of_variation:
 208    BOOST_CHECK_CLOSE(
 209     coefficient_of_variation(dist)
 210     , standard_deviation(dist) / mean(dist), tolerance);
 211    // mode:
 212    BOOST_CHECK_CLOSE(
 213     mode(dist)
 214     , static_cast<RealType>(0.36787944117144232159552377016146L), tolerance);
 215
 216    BOOST_CHECK_CLOSE(
 217     median(dist)
 218     , static_cast<RealType>(exp(dist.location())), tolerance);
 219
 220    BOOST_CHECK_CLOSE(
 221     median(dist),
 222     quantile(dist, static_cast<RealType>(0.5)), tolerance);
 223
 224    // skewness:
 225    BOOST_CHECK_CLOSE(
 226     skewness(dist)
 227     , static_cast<RealType>(729551.38304660255658441529235697L), tolerance);
 228    // kertosis:
 229    BOOST_CHECK_CLOSE(
 230     kurtosis(dist)
 231     , static_cast<RealType>(4312295840576303.2363383232038251L), tolerance);
 232    // kertosis excess:
 233    BOOST_CHECK_CLOSE(
 234     kurtosis_excess(dist)
 235     , static_cast<RealType>(4312295840576300.2363383232038251L), tolerance);
 236
 237    BOOST_CHECK_CLOSE(
 238     range(dist).first
 239     , static_cast<RealType>(0), tolerance);
 240
 241    //
 242    // Special cases:
 243    //
 244    BOOST_CHECK(pdf(dist, 0) == 0);
 245    BOOST_CHECK(cdf(dist, 0) == 0);
 246    BOOST_CHECK(cdf(complement(dist, 0)) == 1);
 247    BOOST_CHECK(quantile(dist, 0) == 0);
 248    BOOST_CHECK(quantile(complement(dist, 1)) == 0);
 249
 250    //
 251    // Error checks:
 252    //
 253    BOOST_MATH_CHECK_THROW(lognormal_distribution<RealType>(0, 0), std::domain_error);
 254    BOOST_MATH_CHECK_THROW(lognormal_distribution<RealType>(2, -1), std::domain_error);
 255    BOOST_MATH_CHECK_THROW(pdf(dist, -1), std::domain_error);
 256    BOOST_MATH_CHECK_THROW(cdf(dist, -1), std::domain_error);
 257    BOOST_MATH_CHECK_THROW(cdf(complement(dist, -1)), std::domain_error);
 258    BOOST_MATH_CHECK_THROW(quantile(dist, 1), std::overflow_error);
 259    BOOST_MATH_CHECK_THROW(quantile(complement(dist, 0)), std::overflow_error);
 260    check_out_of_range<lognormal_distribution<RealType> >(1, 2);
 261
 262 } // template <class RealType>void test_spots(RealType)
 263
 264 BOOST_AUTO_TEST_CASE( test_main )
 265 {
 266
 267    // Check that can generate lognormal distribution using the two convenience methods:
 268    boost::math::lognormal myf1(1., 2); // Using typedef
 269    lognormal_distribution<> myf2(1., 2); // Using default RealType double.
 270
 271   // Test range and support using double only,
 272   // because it supports numeric_limits max for a pseudo-infinity.
 273   BOOST_CHECK_EQUAL(range(myf2).first, 0); // range 0 to +infinity
 274   BOOST_CHECK_EQUAL(range(myf2).second, (std::numeric_limits<double>::max)());
 275   BOOST_CHECK_EQUAL(support(myf2).first, 0); // support 0 to + infinity.
 276   BOOST_CHECK_EQUAL(support(myf2).second, (std::numeric_limits<double>::max)());
 277
 278     // Basic sanity-check spot values.
 279    // (Parameter value, arbitrarily zero, only communicates the floating point type).
 280   test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 %
 281   test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 %
 282 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 283   test_spots(0.0L); // Test long double.
 284 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x0582))
 285   test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
 286 #endif
 287 #else
 288    std::cout << "<note>The long double tests have been disabled on this platform "
 289       "either because the long double overloads of the usual math functions are "
 290       "not available at all, or because they are too inaccurate for these tests "
 291       "to pass.</note>" << std::endl;
 292 #endif
 293
 294
 295 } // BOOST_AUTO_TEST_CASE( test_main )
 296
 297 /*
 298 Running 1 test case...
 299 Tolerance for type float is 0.5 %
 300 Tolerance for type float is 5.96046e-005 %
 301 Tolerance for type float is 5.96046e-005 %
 302 Tolerance for type double is 0.5 %
 303 Tolerance for type double is 1.11022e-013 %
 304 Tolerance for type double is 1.11022e-013 %
 305 Tolerance for type long double is 0.5 %
 306 Tolerance for type long double is 1.11022e-013 %
 307 Tolerance for type long double is 1.11022e-013 %
 308 Tolerance for type class boost::math::concepts::real_concept is 0.5 %
 309 Tolerance for type class boost::math::concepts::real_concept is 1.11022e-013 %
 310 Tolerance for type class boost::math::concepts::real_concept is 1.11022e-013 %
 311 *** No errors detected
 312 */
 313
 314