libs/math/test/test_inverse_gaussian.cpp

   1 // Copyright Paul A. Bristow 2010.
   2 // Copyright John Maddock 2010.
   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 #ifdef _MSC_VER
  10 #  pragma warning (disable : 4224) // nonstandard extension used : formal parameter 'type' was previously defined as a type
  11 // in Boost.test and lexical_cast
  12 #  pragma warning (disable : 4310) // cast truncates constant value
  13 #  pragma warning (disable : 4512) // assignment operator could not be generated
  14
  15 #endif
  16
  17 //#include <pch.hpp> // include directory libs/math/src/tr1/ is needed.
  18
  19 #include <boost/math/tools/test.hpp>
  20 #include <boost/math/concepts/real_concept.hpp> // for real_concept
  21 #define BOOST_TEST_MAIN
  22 #include <boost/test/unit_test.hpp> // Boost.Test
  23 #include <boost/test/tools/floating_point_comparison.hpp>
  24
  25 #include <boost/math/distributions/inverse_gaussian.hpp>
  26 using boost::math::inverse_gaussian_distribution;
  27 using boost::math::inverse_gaussian;
  28
  29 #include <boost/math/tools/test.hpp>
  30 #include "test_out_of_range.hpp"
  31
  32 #include <iostream>
  33 #include <iomanip>
  34 using std::cout;
  35 using std::endl;
  36 using std::setprecision;
  37 #include <limits>
  38 using std::numeric_limits;
  39
  40 template <class RealType>
  41 void check_inverse_gaussian(RealType mean, RealType scale, RealType x, RealType p, RealType q, RealType tol)
  42 {
  43  using boost::math::inverse_gaussian_distribution;
  44
  45   BOOST_CHECK_CLOSE_FRACTION(
  46     ::boost::math::cdf(   // Check cdf
  47     inverse_gaussian_distribution<RealType>(mean, scale),      // distribution.
  48     x),    // random variable.
  49     p,     // probability.
  50     tol);   // tolerance.
  51   BOOST_CHECK_CLOSE_FRACTION(
  52     ::boost::math::cdf( // Check cdf complement
  53     complement(
  54     inverse_gaussian_distribution<RealType>(mean, scale),   // distribution.
  55     x)),   // random variable.
  56     q,      // probability complement.
  57     tol);    // %tolerance.
  58   BOOST_CHECK_CLOSE_FRACTION(
  59     ::boost::math::quantile( // Check quantile
  60     inverse_gaussian_distribution<RealType>(mean, scale),    // distribution.
  61     p),   // probability.
  62     x,   // random variable.
  63     tol);   // tolerance.
  64   BOOST_CHECK_CLOSE_FRACTION(
  65     ::boost::math::quantile( // Check quantile complement
  66     complement(
  67     inverse_gaussian_distribution<RealType>(mean, scale),   // distribution.
  68     q)),   // probability complement.
  69     x,     // random variable.
  70     tol);  // tolerance.
  71
  72    inverse_gaussian_distribution<RealType> dist (mean, scale);
  73
  74    if((p < 0.999) && (q < 0.999))
  75    {  // We can only check this if P is not too close to 1,
  76       // so that we can guarantee Q is accurate:
  77       BOOST_CHECK_CLOSE_FRACTION(
  78         cdf(complement(dist, x)), q, tol); // 1 - cdf
  79       BOOST_CHECK_CLOSE_FRACTION(
  80         quantile(dist, p), x, tol); // quantile(cdf) = x
  81       BOOST_CHECK_CLOSE_FRACTION(
  82         quantile(complement(dist, q)), x, tol); // quantile(complement(1 - cdf)) = x
  83    }
  84 }
  85
  86 template <class RealType>
  87 void test_spots(RealType)
  88 {
  89   // Basic sanity checks
  90   RealType tolerance = static_cast<RealType>(1e-4L); //
  91   cout << "Tolerance for type " << typeid(RealType).name()  << " is " << tolerance << endl;
  92
  93   // Check some bad parameters to the distribution,
  94 #ifndef BOOST_NO_EXCEPTIONS
  95   BOOST_MATH_CHECK_THROW(boost::math::inverse_gaussian_distribution<RealType> nbad1(0, 0), std::domain_error); // zero scale
  96   BOOST_MATH_CHECK_THROW(boost::math::inverse_gaussian_distribution<RealType> nbad1(0, -1), std::domain_error); // negative scale
  97 #else
  98   BOOST_MATH_CHECK_THROW(boost::math::inverse_gaussian_distribution<RealType>(0, 0), std::domain_error); // zero scale
  99   BOOST_MATH_CHECK_THROW(boost::math::inverse_gaussian_distribution<RealType>(0, -1), std::domain_error); // negative scale
 100 #endif
 101
 102   inverse_gaussian_distribution<RealType> w11;
 103
 104   // Error tests:
 105   check_out_of_range<inverse_gaussian_distribution<RealType> >(0.25, 1);
 106
 107   // Check complements.
 108
 109     BOOST_CHECK_CLOSE_FRACTION(
 110      cdf(complement(w11, 1.)), static_cast<RealType>(1) - cdf(w11, 1.), tolerance); // cdf complement
 111     // cdf(complement = 1 - cdf  - but if cdf near unity, then loss of accuracy in cdf,
 112     // but cdf complement is near zero but more accurate.
 113
 114      BOOST_CHECK_CLOSE_FRACTION( // quantile(complement p) == quantile(1 - p)
 115      quantile(complement(w11, static_cast<RealType>(0.5))),
 116      quantile(w11, 1 - static_cast<RealType>(0.5)),
 117      tolerance); // cdf complement
 118
 119   check_inverse_gaussian(
 120      static_cast<RealType>(2),
 121      static_cast<RealType>(3),
 122      static_cast<RealType>(1),
 123      static_cast<RealType>(0.28738674440477374),
 124      static_cast<RealType>(1 - 0.28738674440477374),
 125      tolerance);
 126
 127   RealType tolfeweps = boost::math::tools::epsilon<RealType>() * 5;
 128
 129   inverse_gaussian_distribution<RealType> dist(2, 3);
 130
 131   using namespace std; // ADL of std names.
 132   // mean:
 133   BOOST_CHECK_CLOSE_FRACTION(mean(dist),
 134     static_cast<RealType>(2), tolfeweps);
 135   BOOST_CHECK_CLOSE_FRACTION(scale(dist),
 136     static_cast<RealType>(3), tolfeweps);
 137
 138   // variance:
 139   BOOST_CHECK_CLOSE_FRACTION(variance(dist),
 140     static_cast<RealType>(2.6666666666666666666666666666666666666666666666666666666667L), 1000*tolfeweps);
 141   // std deviation:
 142   BOOST_CHECK_CLOSE_FRACTION(standard_deviation(dist),
 143     static_cast<RealType>(1.632993L), 1000 * tolerance);
 144   //// hazard:
 145   //BOOST_CHECK_CLOSE_FRACTION(hazard(dist, x),
 146   //  pdf(dist, x) / cdf(complement(dist, x)), tolerance);
 147   //// cumulative hazard:
 148   //BOOST_CHECK_CLOSE_FRACTION(chf(dist, x),
 149   //  -log(cdf(complement(dist, x))), tolerance);
 150   // coefficient_of_variation:
 151   BOOST_CHECK_CLOSE_FRACTION(coefficient_of_variation(dist),
 152     standard_deviation(dist) / mean(dist), tolerance);
 153   // mode:
 154   BOOST_CHECK_CLOSE_FRACTION(mode(dist),
 155     static_cast<RealType>(0.8284271L), tolerance);
 156
 157   // median
 158   BOOST_CHECK_CLOSE_FRACTION(median(dist),
 159     static_cast<RealType>(1.5122506636053668L), tolerance);
 160   // Fails for real_concept - because std::numeric_limits<RealType>::digits = 0
 161
 162   // skewness:
 163   BOOST_CHECK_CLOSE_FRACTION(skewness(dist),
 164     static_cast<RealType>(2.449490L), tolerance);
 165   // kurtosis:
 166   BOOST_CHECK_CLOSE_FRACTION(kurtosis(dist),
 167     static_cast<RealType>(10-3), tolerance);
 168   BOOST_CHECK_CLOSE_FRACTION(kurtosis_excess(dist),
 169     static_cast<RealType>(10), tolerance);
 170 } // template <class RealType>void test_spots(RealType)
 171
 172 BOOST_AUTO_TEST_CASE( test_main )
 173 {
 174   using boost::math::inverse_gaussian;
 175   using boost::math::inverse_gaussian_distribution;
 176
 177   //int precision = 17; // std::numeric_limits<double::max_digits10;
 178   double tolfeweps = numeric_limits<double>::epsilon() * 5;
 179   //double tol6decdigits = numeric_limits<float>::epsilon() * 2;
 180   // Check that can generate inverse_gaussian distribution using the two convenience methods:
 181   boost::math::inverse_gaussian w12(1., 2); // Using typedef
 182   inverse_gaussian_distribution<> w23(2., 3); // Using default RealType double.
 183   boost::math::inverse_gaussian w11; // Use default unity values for mean and scale.
 184   // Note NOT myn01() as the compiler will interpret as a function!
 185   BOOST_CHECK_EQUAL(w11.mean(), 1);
 186   BOOST_CHECK_EQUAL(w11.scale(), 1);
 187   BOOST_CHECK_EQUAL(w23.mean(), 2);
 188   BOOST_CHECK_EQUAL(w23.scale(), 3);
 189   BOOST_CHECK_EQUAL(w23.shape(), 1.5L);
 190
 191   // Check the synonyms, provided to allow generic use of find_location and find_scale.
 192   BOOST_CHECK_EQUAL(w11.mean(), w11.location());
 193   BOOST_CHECK_EQUAL(w11.scale(), w11.scale());
 194
 195   BOOST_CHECK_CLOSE_FRACTION(mean(w11), static_cast<double>(1), tolfeweps); // Default mean == unity
 196   BOOST_CHECK_CLOSE_FRACTION(scale(w11), static_cast<double>(1), tolfeweps); // Default mean == unity
 197
 198   // median
 199   // (test double because fails for real_concept because numeric_limits<real_concept>::digits = 0)
 200   BOOST_CHECK_CLOSE_FRACTION(median(w11),
 201     static_cast<double>(0.67584130569523893), tolfeweps);
 202   BOOST_CHECK_CLOSE_FRACTION(median(w23),
 203     static_cast<double>(1.5122506636053668), tolfeweps);
 204
 205   // Initial spot tests using double values from R.
 206   // library(SuppDists)
 207   // formatC(SuppDists::dinverse_gaussian(1, 1, 1), digits=17) ...
 208   BOOST_CHECK_CLOSE_FRACTION( //  x = 1
 209     pdf(w11, 1.), static_cast<double>(0.3989422804014327), tolfeweps); // pdf
 210   BOOST_CHECK_CLOSE_FRACTION(
 211     cdf(w11, 1.), static_cast<double>(0.66810200122317065), 10 * tolfeweps); // cdf
 212
 213   BOOST_CHECK_CLOSE_FRACTION(
 214     pdf(w11, 0.1), static_cast<double>(0.21979480031862672), tolfeweps); // pdf
 215   BOOST_CHECK_CLOSE_FRACTION(
 216     cdf(w11, 0.1), static_cast<double>(0.0040761113207110162), 10 * tolfeweps); // cdf
 217
 218   BOOST_CHECK_CLOSE_FRACTION( // small x
 219     pdf(w11, 0.01), static_cast<double>(2.0811768202028392e-19), tolfeweps); // pdf
 220   BOOST_CHECK_CLOSE_FRACTION(
 221     cdf(w11, 0.01), static_cast<double>(4.122313403318778e-23), 10 * tolfeweps); // cdf
 222
 223   BOOST_CHECK_CLOSE_FRACTION( // smaller x
 224     pdf(w11, 0.001), static_cast<double>(2.4420044378793562e-213),  tolfeweps); // pdf
 225   BOOST_CHECK_CLOSE_FRACTION(
 226     cdf(w11, 0.001), static_cast<double>(4.8791443010851493e-219), 1000 * tolfeweps); // cdf
 227   // 4.8791443010859224e-219 versus 4.8791443010851493e-219 so still 14 decimal digits.
 228
 229   BOOST_CHECK_CLOSE_FRACTION(
 230     quantile(w11, 0.66810200122317065), static_cast<double>(1.), 1 * tolfeweps); // cdf
 231   BOOST_CHECK_CLOSE_FRACTION(
 232     quantile(w11, 0.0040761113207110162), static_cast<double>(0.1), 1 * tolfeweps); // cdf
 233   BOOST_CHECK_CLOSE_FRACTION(
 234     quantile(w11, 4.122313403318778e-23), 0.01, 1 * tolfeweps); // quantile
 235   BOOST_CHECK_CLOSE_FRACTION(
 236     quantile(w11, 2.4420044378793562e-213), 0.001, 0.03); // quantile
 237   // quantile 0.001026926242348481 compared to expected 0.001, so much less accurate,
 238   // but better than R that gives up completely!
 239   // R Error in SuppDists::qinverse_gaussian(4.87914430108515e-219, 1, 1) : Infinite value in NewtonRoot()
 240
 241   BOOST_CHECK_CLOSE_FRACTION(
 242     pdf(w11, 0.5), static_cast<double>(0.87878257893544476), tolfeweps); // pdf
 243   BOOST_CHECK_CLOSE_FRACTION(
 244     cdf(w11, 0.5), static_cast<double>(0.3649755481729598), tolfeweps); // cdf
 245
 246   BOOST_CHECK_CLOSE_FRACTION(
 247     pdf(w11, 2), static_cast<double>(0.10984782236693059), tolfeweps); // pdf
 248   BOOST_CHECK_CLOSE_FRACTION(
 249     cdf(w11, 2), static_cast<double>(.88547542598600637), tolfeweps); // cdf
 250
 251   BOOST_CHECK_CLOSE_FRACTION(
 252     pdf(w11, 10), static_cast<double>(0.00021979480031862676), tolfeweps); // pdf
 253   BOOST_CHECK_CLOSE_FRACTION(
 254     cdf(w11, 10), static_cast<double>(0.99964958546279115), tolfeweps); // cdf
 255
 256   BOOST_CHECK_CLOSE_FRACTION(
 257     pdf(w11, 100), static_cast<double>(2.0811768202028246e-25), tolfeweps); // pdf
 258   BOOST_CHECK_CLOSE_FRACTION(
 259     cdf(w11, 100), static_cast<double>(1), tolfeweps); // cdf
 260   BOOST_CHECK_CLOSE_FRACTION(
 261     pdf(w11, 1000), static_cast<double>(2.4420044378793564e-222), 10 * tolfeweps); // pdf
 262   BOOST_CHECK_CLOSE_FRACTION(
 263     cdf(w11, 1000), static_cast<double>(1.), tolfeweps); // cdf
 264
 265   // A few more misc tests, probably not very useful.
 266   BOOST_CHECK_CLOSE_FRACTION(
 267     cdf(w11, 1.), static_cast<double>(0.66810200122317065), tolfeweps); // cdf
 268   BOOST_CHECK_CLOSE_FRACTION(
 269     cdf(w11, 0.1), static_cast<double>(0.0040761113207110162), tolfeweps * 5); // cdf
 270   // 0.0040761113207110162   0.0040761113207110362
 271   BOOST_CHECK_CLOSE_FRACTION(
 272     cdf(w11, 0.2), static_cast<double>(0.063753567519976254), tolfeweps * 5); // cdf
 273   BOOST_CHECK_CLOSE_FRACTION(
 274     cdf(w11, 0.5), static_cast<double>(0.3649755481729598), tolfeweps); // cdf
 275
 276   BOOST_CHECK_CLOSE_FRACTION(
 277     cdf(w11, 0.9), static_cast<double>(0.62502320258649202), tolfeweps); // cdf
 278   BOOST_CHECK_CLOSE_FRACTION(
 279     cdf(w11, 0.99), static_cast<double>(0.66408247396139031), tolfeweps); // cdf
 280   BOOST_CHECK_CLOSE_FRACTION(
 281     cdf(w11, 0.999), static_cast<double>(0.66770275955311675), tolfeweps); // cdf
 282   BOOST_CHECK_CLOSE_FRACTION(
 283     cdf(w11, 10.), static_cast<double>(0.99964958546279115), tolfeweps); // cdf
 284   BOOST_CHECK_CLOSE_FRACTION(
 285     cdf(w11, 50.), static_cast<double>(0.99999999999992029), tolfeweps); // cdf
 286
 287   BOOST_CHECK_CLOSE_FRACTION(
 288     quantile(w11, 0.3649755481729598), static_cast<double>(0.5), tolfeweps); // quantile
 289   BOOST_CHECK_CLOSE_FRACTION(
 290     quantile(w11, 0.62502320258649202), static_cast<double>(0.9), tolfeweps); // quantile
 291   BOOST_CHECK_CLOSE_FRACTION(
 292     quantile(w11, 0.0040761113207110162), static_cast<double>(0.1), tolfeweps); // quantile
 293
 294   // Wald(2,3) tests
 295   // ===================
 296   BOOST_CHECK_CLOSE_FRACTION( // formatC(SuppDists::dinvGauss(1, 2, 3), digits=17) "0.47490884963330904"
 297     pdf(w23, 1.), static_cast<double>(0.47490884963330904), tolfeweps ); // pdf
 298
 299   BOOST_CHECK_CLOSE_FRACTION(
 300     pdf(w23, 0.1), static_cast<double>(2.8854207087665401e-05), tolfeweps * 2); // pdf
 301   //2.8854207087665452e-005 2.8854207087665401e-005
 302   BOOST_CHECK_CLOSE_FRACTION(
 303     pdf(w23, 10.), static_cast<double>(0.0019822751498574636), tolfeweps); // pdf
 304   BOOST_CHECK_CLOSE_FRACTION(
 305     pdf(w23, 10.), static_cast<double>(0.0019822751498574636), tolfeweps); // pdf
 306
 307   // Bigger changes in mean and scale.
 308
 309   inverse_gaussian w012(0.1, 2);
 310   BOOST_CHECK_CLOSE_FRACTION(
 311     pdf(w012, 1.), static_cast<double>(3.7460367141230404e-36), tolfeweps ); // pdf
 312   BOOST_CHECK_CLOSE_FRACTION(
 313     cdf(w012, 1.), static_cast<double>(1), tolfeweps ); // pdf
 314
 315   inverse_gaussian w0110(0.1, 10);
 316   BOOST_CHECK_CLOSE_FRACTION(
 317     pdf(w0110, 1.), static_cast<double>(1.6279643678071011e-176), 100 * tolfeweps ); // pdf
 318   BOOST_CHECK_CLOSE_FRACTION(
 319     cdf(w0110, 1.), static_cast<double>(1), tolfeweps ); // cdf
 320   BOOST_CHECK_CLOSE_FRACTION(
 321      cdf(complement(w0110, 1.)), static_cast<double>(3.2787685715328683e-179), 1e6 * tolfeweps ); // cdf complement
 322   // Differs because of loss of accuracy.
 323
 324   BOOST_CHECK_CLOSE_FRACTION(
 325     pdf(w0110, 0.1), static_cast<double>(39.894228040143268), tolfeweps ); // pdf
 326   BOOST_CHECK_CLOSE_FRACTION(
 327     cdf(w0110, 0.1), static_cast<double>(0.51989761564832704), 10 * tolfeweps ); // cdf
 328
 329     // Basic sanity-check spot values for all floating-point types..
 330   // (Parameter value, arbitrarily zero, only communicates the floating point type).
 331   test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 %
 332   test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 %
 333 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 334   test_spots(0.0L); // Test long double.
 335 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
 336   test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
 337 #endif
 338 #else
 339   std::cout << "<note>The long double tests have been disabled on this platform "
 340     "either because the long double overloads of the usual math functions are "
 341     "not available at all, or because they are too inaccurate for these tests "
 342     "to pass.</note>" << std::endl;
 343 #endif
 344   /*      */
 345
 346 } // BOOST_AUTO_TEST_CASE( test_main )
 347
 348 /*
 349
 350 Output:
 351
 352
 353 */
 354
 355