ceph/src/boost/libs/math/test/test_geometric.cpp

   1 // test_geometric.cpp
   2
   3 // Copyright Paul A. Bristow 2010.
   4 // Copyright John Maddock 2010.
   5
   6 // Use, modification and distribution are subject to the
   7 // Boost Software License, Version 1.0.
   8 // (See accompanying file LICENSE_1_0.txt
   9 // or copy at http://www.boost.org/LICENSE_1_0.txt)
  10
  11 // Tests for Geometric Distribution.
  12
  13 // Note that these defines must be placed BEFORE #includes.
  14 #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
  15 // because several tests overflow & underflow by design.
  16 #define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
  17
  18 #ifdef _MSC_VER
  19 #  pragma warning(disable: 4127) // conditional expression is constant.
  20 #endif
  21
  22 #if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
  23 #  define TEST_FLOAT
  24 #  define TEST_DOUBLE
  25 #  define TEST_LDOUBLE
  26 #  define TEST_REAL_CONCEPT
  27 #endif
  28
  29 #include <boost/math/tools/test.hpp>
  30 #include <boost/math/concepts/real_concept.hpp> // for real_concept
  31 using ::boost::math::concepts::real_concept;
  32
  33 #include <boost/math/distributions/geometric.hpp> // for geometric_distribution
  34 using boost::math::geometric_distribution;
  35 using boost::math::geometric; // using typedef for geometric_distribution<double>
  36
  37 #include <boost/math/distributions/negative_binomial.hpp> // for some comparisons.
  38
  39 #define BOOST_TEST_MAIN
  40 #include <boost/test/unit_test.hpp> // for test_main
  41 #include <boost/test/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE_FRACTION
  42 #include "test_out_of_range.hpp"
  43
  44 #include <iostream>
  45 using std::cout;
  46 using std::endl;
  47 using std::setprecision;
  48 using std::showpoint;
  49 #include <limits>
  50 using std::numeric_limits;
  51
  52 template <class RealType>
  53 void test_spot( // Test a single spot value against 'known good' values.
  54                RealType k,    // Number of failures.
  55                RealType p,    // Probability of success_fraction.
  56                RealType P,    // CDF probability.
  57                RealType Q,    // Complement of CDF.
  58                RealType tol)  // Test tolerance.
  59 {
  60    boost::math::geometric_distribution<RealType> g(p);
  61    BOOST_CHECK_EQUAL(p, g.success_fraction());
  62    BOOST_CHECK_CLOSE_FRACTION(cdf(g, k), P, tol);
  63
  64   if((P < 0.99) && (Q < 0.99))
  65   {
  66     // We can only check this if P is not too close to 1,
  67     // so that we can guarantee that Q is free of error:
  68     //
  69     BOOST_CHECK_CLOSE_FRACTION(
  70       cdf(complement(g, k)), Q, tol);
  71     if(k != 0)
  72     {
  73       BOOST_CHECK_CLOSE_FRACTION(
  74         quantile(g, P), k, tol);
  75     }
  76     else
  77     {
  78       // Just check quantile is very small:
  79       if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent)
  80         && (boost::is_floating_point<RealType>::value))
  81       {
  82         // Limit where this is checked: if exponent range is very large we may
  83         // run out of iterations in our root finding algorithm.
  84         BOOST_CHECK(quantile(g, P) < boost::math::tools::epsilon<RealType>() * 10);
  85       }
  86     }
  87     if(k != 0)
  88     {
  89       BOOST_CHECK_CLOSE_FRACTION(
  90         quantile(complement(g, Q)), k, tol);
  91     }
  92     else
  93     {
  94       // Just check quantile is very small:
  95       if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent)
  96         && (boost::is_floating_point<RealType>::value))
  97       {
  98         // Limit where this is checked: if exponent range is very large we may
  99         // run out of iterations in our root finding algorithm.
 100         BOOST_CHECK(quantile(complement(g, Q)) < boost::math::tools::epsilon<RealType>() * 10);
 101       }
 102     }
 103   } //   if((P < 0.99) && (Q < 0.99))
 104
 105     // Parameter estimation test:  estimate success ratio:
 106     BOOST_CHECK_CLOSE_FRACTION(
 107       geometric_distribution<RealType>::find_lower_bound_on_p(
 108       1+k, P),
 109       p, 0.02); // Wide tolerance needed for some tests.
 110    // Note we bump up the sample size here, purely for the sake of the test,
 111     // internally the function has to adjust the sample size so that we get
 112     // the right upper bound, our test undoes this, so we can verify the result.
 113     BOOST_CHECK_CLOSE_FRACTION(
 114       geometric_distribution<RealType>::find_upper_bound_on_p(
 115       1+k+1, Q),
 116       p, 0.02);
 117
 118     if(Q < P)
 119     {
 120        //
 121        // We check two things here, that the upper and lower bounds
 122        // are the right way around, and that they do actually bracket
 123        // the naive estimate of p = successes / (sample size)
 124        //
 125       BOOST_CHECK(
 126         geometric_distribution<RealType>::find_lower_bound_on_p(
 127         1+k, Q)
 128         <=
 129         geometric_distribution<RealType>::find_upper_bound_on_p(
 130         1+k, Q)
 131         );
 132       BOOST_CHECK(
 133         geometric_distribution<RealType>::find_lower_bound_on_p(
 134         1+k, Q)
 135         <=
 136         1 / (1+k)
 137         );
 138       BOOST_CHECK(
 139         1 / (1+k)
 140         <=
 141         geometric_distribution<RealType>::find_upper_bound_on_p(
 142         1+k, Q)
 143         );
 144     }
 145     else
 146     {
 147        // As above but when P is small.
 148       BOOST_CHECK(
 149         geometric_distribution<RealType>::find_lower_bound_on_p(
 150         1+k, P)
 151         <=
 152         geometric_distribution<RealType>::find_upper_bound_on_p(
 153         1+k, P)
 154         );
 155       BOOST_CHECK(
 156         geometric_distribution<RealType>::find_lower_bound_on_p(
 157         1+k,  P)
 158         <=
 159         1 / (1+k)
 160         );
 161       BOOST_CHECK(
 162         1 / (1+k)
 163         <=
 164         geometric_distribution<RealType>::find_upper_bound_on_p(
 165         1+k, P)
 166         );
 167     }
 168
 169     // Estimate sample size:
 170     BOOST_CHECK_CLOSE_FRACTION(
 171       geometric_distribution<RealType>::find_minimum_number_of_trials(
 172       k, p, P),
 173       1+k, 0.02); // Can differ 50 to 51 for small p
 174     BOOST_CHECK_CLOSE_FRACTION(
 175       geometric_distribution<RealType>::find_maximum_number_of_trials(
 176          k, p, Q),
 177       1+k, 0.02);
 178
 179 } // test_spot
 180
 181 template <class RealType> // Any floating-point type RealType.
 182 void test_spots(RealType)
 183 {
 184   // Basic sanity checks.
 185   // Most test data is to double precision (17 decimal digits) only,
 186
 187   cout << "Floating point Type is " << typeid(RealType).name() << endl;
 188
 189   // so set tolerance to 1000 eps expressed as a fraction,
 190   // or 1000 eps of type double expressed as a fraction,
 191   // whichever is the larger.
 192
 193   RealType tolerance = (std::max)
 194     (boost::math::tools::epsilon<RealType>(),
 195     static_cast<RealType>(std::numeric_limits<double>::epsilon()));
 196   tolerance *= 10; // 10 eps
 197
 198   cout << "Tolerance = " << tolerance << "." << endl;
 199
 200   RealType tol1eps = boost::math::tools::epsilon<RealType>(); // Very tight, suit exact values.
 201   //RealType tol2eps = boost::math::tools::epsilon<RealType>() * 2; // Tight,  values.
 202   RealType tol5eps = boost::math::tools::epsilon<RealType>() * 5; // Wider 5 epsilon.
 203   cout << "Tolerance 5 eps = " << tol5eps << "." << endl;
 204
 205
 206   // Sources of spot test values are mainly R.
 207
 208   using boost::math::geometric_distribution;
 209   using boost::math::geometric;
 210   using boost::math::cdf;
 211   using boost::math::pdf;
 212   using boost::math::quantile;
 213   using boost::math::complement;
 214
 215   BOOST_MATH_STD_USING // for std math functions
 216
 217   // Test geometric using cdf spot values R
 218   // These test quantiles and complements as well.
 219
 220   test_spot(  //
 221   static_cast<RealType>(2),   // Number of failures, k
 222   static_cast<RealType>(0.5), // Probability of success as fraction, p
 223   static_cast<RealType>(0.875L), // Probability of result (CDF), P
 224   static_cast<RealType>(0.125L),  // complement CCDF Q = 1 - P
 225   tolerance);
 226
 227   test_spot( //
 228   static_cast<RealType>(0),    // Number of failures, k
 229   static_cast<RealType>(0.25), // Probability of success as fraction, p
 230   static_cast<RealType>(0.25),   // Probability of result (CDF), P
 231   static_cast<RealType>(0.75),   // Q = 1 - P
 232   tolerance);
 233
 234   test_spot(
 235     // R formatC(pgeom(10,0.25), digits=17) [1] "0.95776486396789551"
 236     // formatC(pgeom(10,0.25, FALSE), digits=17) [1] "0.042235136032104499"
 237
 238   static_cast<RealType>(10),  // Number of failures, k
 239   static_cast<RealType>(0.25),  // Probability of success, p
 240   static_cast<RealType>(0.95776486396789551L),  // Probability of result (CDF), P
 241   static_cast<RealType>(0.042235136032104499L), // Q = 1 - P
 242   tolerance);
 243
 244   test_spot(  //
 245   // > R formatC(pgeom(50,0.25, TRUE), digits=17) [1] "0.99999957525875771"
 246   // > R formatC(pgeom(50,0.25, FALSE), digits=17) [1] "4.2474124232020353e-07"
 247   static_cast<RealType>(50),     // Number of failures, k
 248   static_cast<RealType>(0.25),     // Probability of success, p
 249   static_cast<RealType>(0.99999957525875771),  // Probability of result (CDF), P
 250   static_cast<RealType>(4.2474124232020353e-07),   // Q = 1 - P
 251   tolerance);
 252   /*
 253   // This causes failures in find_upper_bound_on_p p is small branch.
 254   test_spot(  // formatC(pgeom(50,0.01, TRUE), digits=17)[1] "0.40104399353383874"
 255     // > formatC(pgeom(50,0.01, FALSE), digits=17) [1] "0.59895600646616121"
 256   static_cast<RealType>(50), // Number of failures, k
 257   static_cast<RealType>(0.01),   // Probability of success, p
 258   static_cast<RealType>(0.40104399353383874),   // Probability of result (CDF), P
 259   static_cast<RealType>(0.59895600646616121),   // Q = 1 - P
 260   tolerance);
 261   */
 262
 263   test_spot( // > formatC(pgeom(50,0.99, TRUE), digits=17) [1] "                 1"
 264     // formatC(pgeom(50,0.99, FALSE), digits=17) [1] "1.0000000000000364e-102"
 265   static_cast<RealType>(50),     // Number of failures, k
 266   static_cast<RealType>(0.99),    // Probability of success, p
 267   static_cast<RealType>(1), // Probability of result (CDF), P
 268   static_cast<RealType>(1.0000000000000364e-102),   // Q = 1 - P
 269   tolerance);
 270
 271   test_spot(  // > formatC(pgeom(1,0.99, TRUE), digits=17) [1] "0.99990000000000001"
 272     // > formatC(pgeom(1,0.99, FALSE), digits=17) [1] "0.00010000000000000009"
 273   static_cast<RealType>(1),     // Number of failures, k
 274   static_cast<RealType>(0.99),                    // Probability of success, p
 275   static_cast<RealType>(0.9999),     // Probability of result (CDF), P
 276   static_cast<RealType>(0.0001),   // Q = 1 - P
 277   tolerance);
 278
 279 if(std::numeric_limits<RealType>::is_specialized)
 280 { // An extreme value test that is more accurate than using negative binomial.
 281   // Since geometric only uses exp and log functions.
 282   test_spot(  // > formatC(pgeom(10000, 0.001, TRUE), digits=17) [1] "0.99995487182736897"
 283 // > formatC(pgeom(10000,0.001, FALSE), digits=17) [1] "4.5128172631071587e-05"
 284   static_cast<RealType>(10000L), // Number of failures, k
 285   static_cast<RealType>(0.001L),                    // Probability of success, p
 286   static_cast<RealType>(0.99995487182736897L),     // Probability of result (CDF), P
 287   static_cast<RealType>(4.5128172631071587e-05L),   // Q = 1 - P
 288   tolerance); //
 289   } // numeric_limit is specialized
 290  // End of single spot tests using RealType
 291
 292   // Tests on PDF:
 293
 294   BOOST_CHECK_CLOSE_FRACTION( //> formatC(dgeom(0,0.5), digits=17)[1] " 0.5"
 295   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
 296   static_cast<RealType>(0.0) ),  // Number of failures, k is very small but not integral,
 297   static_cast<RealType>(0.5), // nearly success probability.
 298   tolerance);
 299
 300   BOOST_CHECK_CLOSE_FRACTION( //> formatC(dgeom(0,0.5), digits=17)[1] "    0.5"
 301     //  R treates geom as a discrete distribution.
 302     // > formatC(dgeom(1.999999,0.5, FALSE), digits=17) [1] "   0"
 303     // Warning message:
 304     // In dgeom(1.999999, 0.5, FALSE) : non-integer x = 1.999999
 305   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
 306   static_cast<RealType>(0.0001L) ),  // Number of failures, k is very small but not integral,
 307   static_cast<RealType>(0.4999653438420768L), // nearly success probability.
 308   tolerance);
 309
 310   BOOST_CHECK_CLOSE_FRACTION( // > formatC(pgeom(0.0001,0.5, TRUE), digits=17)[1] " 0.5"
 311     // > formatC(pgeom(0.0001,0.5, FALSE), digits=17) [1] "               0.5"
 312     //  R treates geom as a discrete distribution.
 313   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
 314   static_cast<RealType>(0.0001L) ),  // Number of failures, k is very small but not integral,
 315   static_cast<RealType>(0.4999653438420768L), // nearly success probability.
 316   tolerance);
 317
 318   BOOST_CHECK_CLOSE_FRACTION( // formatC(dgeom(1,0.01), digits=17)[1] "0.0099000000000000008"
 319   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.01L)),
 320   static_cast<RealType>(1) ),  // Number of failures, k
 321   static_cast<RealType>(0.0099000000000000008), //
 322   tolerance);
 323
 324   BOOST_CHECK_CLOSE_FRACTION( //> formatC(dgeom(1,0.99), digits=17)[1] "0.0099000000000000043"
 325   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.99L)),
 326   static_cast<RealType>(1) ),  // Number of failures, k
 327   static_cast<RealType>(0.00990000000000000043L), //
 328   tolerance);
 329
 330   BOOST_CHECK_CLOSE_FRACTION( //> > formatC(dgeom(0,0.99), digits=17)[1] "0.98999999999999999"
 331   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.99L)),
 332   static_cast<RealType>(0) ),  // Number of failures, k
 333   static_cast<RealType>(0.98999999999999999L), //
 334   tolerance);
 335
 336   // p  near unity.
 337   BOOST_CHECK_CLOSE_FRACTION( // > formatC(dgeom(100,0.99), digits=17)[1] "9.9000000000003448e-201"
 338   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.99L)),
 339   static_cast<RealType>(100) ),  // Number of failures, k
 340   static_cast<RealType>(9.9000000000003448e-201L), //
 341   100 * tolerance); // Note difference
 342
 343     // p nearer unity.
 344   BOOST_CHECK_CLOSE_FRACTION( //
 345   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999)),
 346   static_cast<RealType>(10) ),  // Number of failures, k
 347   // static_cast<double>(9.9989999999889024e-41), // Boost.Math
 348   // static_cast<float>(1.00156406e-040)
 349   static_cast<RealType>(9.999e-41), // exact from 100 digit calculator.
 350   2e3 * tolerance); // Note bigger tolerance needed.
 351
 352   // Moshier Cephes 100 digits calculator says 9.999e-41
 353   //0.9999*pow(1-0.9999,10)
 354   // 9.9990000000000000000000000000000000000000000000000000000000000000000000E-41
 355   // 9.998999999988988e-041
 356   // > formatC(dgeom(10, 0.9999), digits=17) [1] "9.9989999999889024e-41"
 357   // p *  pow(q, k)         9.9989999999889880e-041
 358   // exp(p * k * log1p(-p)) 9.9989999999889024e-041
 359
 360
 361
 362   // 0.9999999999 * pow(1-0.9999999999,10)=  9.9999999990E-101
 363   // > formatC(dgeom(10,0.9999999999), digits=17)  [1] "1.0000008273040127e-100"
 364   BOOST_CHECK_CLOSE_FRACTION( //
 365   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999999999L)),
 366   static_cast<RealType>(10) ),  //
 367   static_cast<RealType>(9.9999999990E-101L), // 1.0000008273040179e-100
 368   1e9 * tolerance); // Note big tolerance needed.
 369   // 1.0000008273040179e-100  Boost.Math
 370   // 1.0000008273040127e-100  R
 371   // 0.9999999990000004e-100  100 digit calculator 'exact'
 372
 373   BOOST_CHECK_CLOSE_FRACTION( //
 374   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.00000000001L)),
 375   static_cast<RealType>(10) ),  //
 376   static_cast<RealType>(9.999999999e-12L), // get 9.9999999989999994e-012
 377   1 * tolerance); // Note small tolerance needed.
 378
 379
 380     BOOST_CHECK_CLOSE_FRACTION( //
 381   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.00000000001L)),
 382   static_cast<RealType>(1000) ),  //
 383   static_cast<RealType>(9.9999999e-12L), // get 9.9999998999999913e-012
 384   tolerance); // Note small tolerance needed.
 385
 386
 387   ///////////////////////////////////////////////////
 388   BOOST_CHECK_CLOSE_FRACTION( //
 389     // > formatC(dgeom(0.0001,0.5, FALSE), digits=17) [1] "               0.5"
 390     //  R treates geom as a discrete distribution.
 391     // But Boost.Math is continuous, so if you want R behaviour,
 392     // make number of failures, k into an integer with the floor function.
 393   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
 394   static_cast<RealType>(floor(0.0001L)) ),  // Number of failures, k is very small but MADE integral,
 395   static_cast<RealType>(0.5), // nearly success probability.
 396   tolerance);
 397
 398   // R switches over at about 1e7 from k = 0, returning 0.5,  to k = 1, returning 0.25.
 399   // Boost.Math does not do this, even for 0.9999999999999999
 400   // > formatC(pgeom(0.999999,0.5, FALSE), digits=17) [1] "               0.5"
 401   // > formatC(pgeom(0.9999999,0.5, FALSE), digits=17) [1] "              0.25"
 402
 403   BOOST_CHECK_CLOSE_FRACTION( // > formatC(pgeom(0.0001,0.5, TRUE), digits=17)[1] "               0.5"
 404     // > formatC(pgeom(0.0001,0.5, FALSE), digits=17) [1] "               0.5"
 405     //  R treates geom as a discrete distribution.
 406     // But Boost.Math is continuous, so if you want R behaviour,
 407     // make number of failures, k into an integer with the floor function.
 408   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
 409   static_cast<RealType>(floor(0.9999999999999999L)) ),  // Number of failures, k is very small but MADE integral,
 410   static_cast<RealType>(0.5), // nearly success probability.
 411   tolerance);
 412
 413   BOOST_CHECK_CLOSE_FRACTION( // > formatC(pgeom(0.0001,0.5, TRUE), digits=17)[1] "               0.5"
 414     // > formatC(pgeom(0.0001,0.5, FALSE), digits=17) [1] "               0.5"
 415     //  R treates geom as a discrete distribution.
 416     // But Boost.Math is continuous, so if you want R behaviour,
 417     // make number of failures, k into an integer with the floor function.
 418   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
 419   static_cast<RealType>(floor(1. - tolerance)) ),
 420   // Number of failures, k is very small but MADE integral,
 421   // Need to use tolerance here,
 422   // as epsilon is ill-defined for Real concept:
 423   // numeric_limits<RealType>::epsilon()  0
 424   static_cast<RealType>(0.5), // nearly success probability.
 425   tolerance * 10);
 426
 427   BOOST_CHECK_CLOSE_FRACTION(
 428   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.0001L)),
 429   static_cast<RealType>(2)),  // k = 2.
 430   static_cast<RealType>(9.99800010e-5L), // 'exact '
 431   tolerance);
 432
 433   //> formatC(dgeom(2, 0.9999), digits=17) [1] "9.9989999999977806e-09"
 434   BOOST_CHECK_CLOSE_FRACTION(
 435   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999L)),
 436   static_cast<RealType>(2)),  // k = 0
 437   static_cast<RealType>(9.999e-9L), // 'exact'
 438   1000*tolerance);
 439
 440   BOOST_CHECK_CLOSE_FRACTION(
 441   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999L)),
 442   static_cast<RealType>(3)),  // k = 3
 443   static_cast<RealType>(9.999e-13L), // get
 444   1000*tolerance);
 445
 446   BOOST_CHECK_CLOSE_FRACTION(
 447   pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999L)),
 448   static_cast<RealType>(5)),  // k = 5
 449   static_cast<RealType>(9.999e-21L), //  9.9989999999944947e-021
 450   1000*tolerance);
 451
 452
 453   BOOST_CHECK_CLOSE_FRACTION(
 454   pdf(geometric_distribution<RealType>( static_cast<RealType>(0.0001L)),
 455   static_cast<RealType>(3)),  // k = 0.
 456   static_cast<RealType>(9.99700029999e-5L), //
 457   tolerance);
 458    // Tests on cdf:
 459   // MathCAD pgeom k, r, p) == failures, successes, probability.
 460
 461   BOOST_CHECK_CLOSE_FRACTION(cdf(
 462     geometric_distribution<RealType>(static_cast<RealType>(0.5)), // prob 0.5
 463     static_cast<RealType>(0) ), // k = 0
 464     static_cast<RealType>(0.5), // probability =p
 465     tolerance);
 466
 467   BOOST_CHECK_CLOSE_FRACTION(cdf(complement(
 468     geometric_distribution<RealType>(static_cast<RealType>(0.5)), //
 469     static_cast<RealType>(0) )), // k = 0
 470     static_cast<RealType>(0.5), // probability =
 471     tolerance);
 472
 473   BOOST_CHECK_CLOSE_FRACTION(cdf(
 474     geometric_distribution<RealType>(static_cast<RealType>(0.25)), // prob 0.5
 475     static_cast<RealType>(1) ), // k = 0
 476     static_cast<RealType>(0.4375L), // probability =p
 477     tolerance);
 478
 479   BOOST_CHECK_CLOSE_FRACTION(cdf(complement(
 480     geometric_distribution<RealType>(static_cast<RealType>(0.25)), //
 481     static_cast<RealType>(1) )), // k = 0
 482     static_cast<RealType>(1-0.4375L), // probability =
 483     tolerance);
 484
 485   BOOST_CHECK_CLOSE_FRACTION(cdf(complement(
 486     geometric_distribution<RealType>(static_cast<RealType>(0.5)), //
 487     static_cast<RealType>(1) )), // k = 0
 488     static_cast<RealType>(0.25), // probability = exact 0.25
 489     tolerance);
 490
 491   BOOST_CHECK_CLOSE_FRACTION( //
 492     cdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
 493     static_cast<RealType>(4)),  // k =4.
 494     static_cast<RealType>(0.96875L), // exact
 495     tolerance);
 496
 497
 498   // Tests of other functions, mean and other moments ...
 499
 500   geometric_distribution<RealType> dist(static_cast<RealType>(0.25));
 501   // mean:
 502   BOOST_CHECK_CLOSE_FRACTION(
 503     mean(dist), static_cast<RealType>((1 - 0.25) /0.25), tol5eps);
 504   BOOST_CHECK_CLOSE_FRACTION(
 505     mode(dist), static_cast<RealType>(0), tol1eps);
 506   // variance:
 507   BOOST_CHECK_CLOSE_FRACTION(
 508     variance(dist), static_cast<RealType>((1 - 0.25) / (0.25 * 0.25)), tol5eps);
 509
 510   // std deviation:
 511   // sqrt(0.75/0.125)
 512
 513   BOOST_CHECK_CLOSE_FRACTION(
 514     standard_deviation(dist), //
 515     static_cast<RealType>(sqrt((1.0L - 0.25L) / (0.25L * 0.25L))), // using 100 digit calc
 516     tol5eps);
 517
 518   BOOST_CHECK_CLOSE_FRACTION(
 519     skewness(dist), //
 520     static_cast<RealType>((2-0.25L) /sqrt(0.75L)),
 521     // using calculator
 522     tol5eps);
 523   BOOST_CHECK_CLOSE_FRACTION(
 524     kurtosis_excess(dist), //
 525     static_cast<RealType>(6 + 0.0625L/0.75L), //
 526     tol5eps);
 527   // 6.083333333333333  6.166666666666667
 528   BOOST_CHECK_CLOSE_FRACTION(
 529     kurtosis(dist), // true
 530     static_cast<RealType>(9 + 0.0625L/0.75L), //
 531     tol5eps);
 532   // hazard:
 533   RealType x = static_cast<RealType>(0.125);
 534   BOOST_CHECK_CLOSE_FRACTION(
 535   hazard(dist, x)
 536   , pdf(dist, x) / cdf(complement(dist, x)), tol5eps);
 537   // cumulative hazard:
 538   BOOST_CHECK_CLOSE_FRACTION(
 539   chf(dist, x), -log(cdf(complement(dist, x))), tol5eps);
 540   // coefficient_of_variation:
 541   BOOST_CHECK_CLOSE_FRACTION(
 542   coefficient_of_variation(dist)
 543   , standard_deviation(dist) / mean(dist), tol5eps);
 544
 545   // Special cases for PDF:
 546   BOOST_CHECK_EQUAL(
 547   pdf(
 548   geometric_distribution<RealType>(static_cast<RealType>(0)), //
 549   static_cast<RealType>(0)),
 550   static_cast<RealType>(0) );
 551
 552   BOOST_CHECK_EQUAL(
 553   pdf(
 554   geometric_distribution<RealType>(static_cast<RealType>(0)),
 555   static_cast<RealType>(0.0001)),
 556   static_cast<RealType>(0) );
 557
 558   BOOST_CHECK_EQUAL(
 559   pdf(
 560   geometric_distribution<RealType>(static_cast<RealType>(1)),
 561   static_cast<RealType>(0.001)),
 562   static_cast<RealType>(0) );
 563
 564   BOOST_CHECK_EQUAL(
 565   pdf(
 566   geometric_distribution<RealType>(static_cast<RealType>(1)),
 567   static_cast<RealType>(8)),
 568   static_cast<RealType>(0) );
 569
 570   BOOST_CHECK_SMALL(
 571   pdf(
 572    geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 573   static_cast<RealType>(0))-
 574   static_cast<RealType>(0.25),
 575   2 * boost::math::tools::epsilon<RealType>() ); // Expect exact, but not quite.
 576   // numeric_limits<RealType>::epsilon()); // Not suitable for real concept!
 577
 578   // Quantile boundary cases checks:
 579   BOOST_CHECK_EQUAL(
 580   quantile(  // zero P < cdf(0) so should be exactly zero.
 581   geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 582   static_cast<RealType>(0)),
 583   static_cast<RealType>(0));
 584
 585   BOOST_CHECK_EQUAL(
 586   quantile(  // min P < cdf(0) so should be exactly zero.
 587   geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 588   static_cast<RealType>(boost::math::tools::min_value<RealType>())),
 589   static_cast<RealType>(0));
 590
 591   BOOST_CHECK_CLOSE_FRACTION(
 592   quantile(  // Small P < cdf(0) so should be near zero.
 593   geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 594   static_cast<RealType>(boost::math::tools::epsilon<RealType>())), //
 595   static_cast<RealType>(0),
 596     tol5eps);
 597
 598   BOOST_CHECK_CLOSE_FRACTION(
 599   quantile(  // Small P < cdf(0) so should be exactly zero.
 600   geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 601   static_cast<RealType>(0.0001)),
 602   static_cast<RealType>(0),
 603     tolerance);
 604
 605   //BOOST_CHECK(  // Fails with overflow for real_concept
 606   //quantile(  // Small P near 1 so k failures should be big.
 607   //geometric_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
 608   //static_cast<RealType>(1 - boost::math::tools::epsilon<RealType>())) <=
 609   //static_cast<RealType>(189.56999032670058)  // 106.462769 for float
 610   //);
 611
 612   if(std::numeric_limits<RealType>::has_infinity)
 613   { // BOOST_CHECK tests for infinity using std::numeric_limits<>::infinity()
 614     // Note that infinity is not implemented for real_concept, so these tests
 615     // are only done for types, like built-in float, double.. that have infinity.
 616     // Note that these assume that  BOOST_MATH_OVERFLOW_ERROR_POLICY is NOT throw_on_error.
 617     // #define BOOST_MATH_THROW_ON_OVERFLOW_POLICY ==  throw_on_error would throw here.
 618     // #define BOOST_MAT_DOMAIN_ERROR_POLICY IS defined throw_on_error,
 619     //  so the throw path of error handling is tested below with BOOST_MATH_CHECK_THROW tests.
 620
 621     BOOST_CHECK(
 622     quantile(  // At P == 1 so k failures should be infinite.
 623     geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 624     static_cast<RealType>(1)) ==
 625     //static_cast<RealType>(boost::math::tools::infinity<RealType>())
 626     static_cast<RealType>(std::numeric_limits<RealType>::infinity()) );
 627
 628     BOOST_CHECK_EQUAL(
 629     quantile(  // At 1 == P  so should be infinite.
 630     geometric_distribution<RealType>( static_cast<RealType>(0.25)),
 631     static_cast<RealType>(1)), //
 632     std::numeric_limits<RealType>::infinity() );
 633
 634     BOOST_CHECK_EQUAL(
 635     quantile(complement(  // Q zero 1 so P == 1 < cdf(0) so should be exactly infinity.
 636     geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 637     static_cast<RealType>(0))),
 638     std::numeric_limits<RealType>::infinity() );
 639    } // test for infinity using std::numeric_limits<>::infinity()
 640   else
 641   { // real_concept case, so check it throws rather than returning infinity.
 642     BOOST_CHECK_EQUAL(
 643     quantile(  // At P == 1 so k failures should be infinite.
 644     geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 645     static_cast<RealType>(1)),
 646     boost::math::tools::max_value<RealType>() );
 647
 648     BOOST_CHECK_EQUAL(
 649     quantile(complement(  // Q zero 1 so P == 1 < cdf(0) so should be exactly infinity.
 650     geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 651     static_cast<RealType>(0))),
 652     boost::math::tools::max_value<RealType>());
 653   } // has infinity
 654
 655   BOOST_CHECK( // Should work for built-in and real_concept.
 656   quantile(complement(  // Q near to 1 so P nearly 1, so should be large > 300.
 657   geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 658   static_cast<RealType>(boost::math::tools::min_value<RealType>())))
 659    >= static_cast<RealType>(300) );
 660
 661   BOOST_CHECK_EQUAL(
 662   quantile(  //  P ==  0 < cdf(0) so should be zero.
 663   geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 664   static_cast<RealType>(0)),
 665   static_cast<RealType>(0));
 666
 667   // Quantile Complement boundary cases:
 668
 669   BOOST_CHECK_EQUAL(
 670   quantile(complement(  // Q = 1 so P = 0 < cdf(0) so should be exactly zero.
 671   geometric_distribution<RealType>( static_cast<RealType>(0.25)),
 672   static_cast<RealType>(1))),
 673   static_cast<RealType>(0)
 674   );
 675
 676   BOOST_CHECK_EQUAL(
 677   quantile(complement(  // Q very near 1 so P == epsilon < cdf(0) so should be exactly zero.
 678   geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 679   static_cast<RealType>(1 - boost::math::tools::epsilon<RealType>()))),
 680   static_cast<RealType>(0)
 681   );
 682
 683   // Check that duff arguments throw domain_error:
 684
 685   BOOST_MATH_CHECK_THROW(
 686   pdf( // Negative success_fraction!
 687   geometric_distribution<RealType>(static_cast<RealType>(-0.25)),
 688   static_cast<RealType>(0)), std::domain_error);
 689   BOOST_MATH_CHECK_THROW(
 690   pdf( // Success_fraction > 1!
 691   geometric_distribution<RealType>(static_cast<RealType>(1.25)),
 692   static_cast<RealType>(0)),
 693   std::domain_error);
 694   BOOST_MATH_CHECK_THROW(
 695   pdf( // Negative k argument !
 696   geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 697   static_cast<RealType>(-1)),
 698   std::domain_error);
 699   //BOOST_MATH_CHECK_THROW(
 700   //pdf( // check limit on k (failures)
 701   //geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 702   //std::numeric_limits<RealType>infinity()),
 703   //std::domain_error);
 704   BOOST_MATH_CHECK_THROW(
 705   cdf(  // Negative k argument !
 706   geometric_distribution<RealType>(static_cast<RealType>(0.25)),
 707   static_cast<RealType>(-1)),
 708   std::domain_error);
 709   BOOST_MATH_CHECK_THROW(
 710   cdf( // Negative success_fraction!
 711   geometric_distribution<RealType>(static_cast<RealType>(-0.25)),
 712   static_cast<RealType>(0)), std::domain_error);
 713   BOOST_MATH_CHECK_THROW(
 714   cdf( // Success_fraction > 1!
 715   geometric_distribution<RealType>(static_cast<RealType>(1.25)),
 716   static_cast<RealType>(0)), std::domain_error);
 717   BOOST_MATH_CHECK_THROW(
 718   quantile(  // Negative success_fraction!
 719   geometric_distribution<RealType>(static_cast<RealType>(-0.25)),
 720   static_cast<RealType>(0)), std::domain_error);
 721   BOOST_MATH_CHECK_THROW(
 722   quantile( // Success_fraction > 1!
 723   geometric_distribution<RealType>(static_cast<RealType>(1.25)),
 724   static_cast<RealType>(0)), std::domain_error);
 725    check_out_of_range<geometric_distribution<RealType> >(0.5);
 726   // End of check throwing 'duff' out-of-domain values.
 727
 728   { // Compare geometric and negative binomial functions.
 729     using boost::math::negative_binomial_distribution;
 730     using boost::math::geometric_distribution;
 731
 732     RealType k = static_cast<RealType>(2.L);
 733     RealType alpha = static_cast<RealType>(0.05L);
 734     RealType p = static_cast<RealType>(0.5L);
 735
 736     BOOST_CHECK_CLOSE_FRACTION( // Successes parameter in negative binomial is 1 for geometric.
 737       geometric_distribution<RealType>::find_lower_bound_on_p(k, alpha),
 738       negative_binomial_distribution<RealType>::find_lower_bound_on_p(k, static_cast<RealType>(1), alpha),
 739       tolerance);
 740     BOOST_CHECK_CLOSE_FRACTION( // Successes parameter in negative binomial is 1 for geometric.
 741       geometric_distribution<RealType>::find_upper_bound_on_p(k, alpha),
 742       negative_binomial_distribution<RealType>::find_upper_bound_on_p(k, static_cast<RealType>(1), alpha),
 743       tolerance);
 744     BOOST_CHECK_CLOSE_FRACTION( // Should be identical - successes parameter is not used.
 745        geometric_distribution<RealType>::find_maximum_number_of_trials(k, p, alpha),
 746       negative_binomial_distribution<RealType>::find_maximum_number_of_trials(k, p, alpha),
 747     tolerance);
 748   }
 749     //geometric::find_upper_bound_on_p(k, alpha);
 750    return;
 751 } // template <class RealType> void test_spots(RealType) // Any floating-point type RealType.
 752
 753 BOOST_AUTO_TEST_CASE( test_main )
 754 {
 755   // Check that can generate geometric distribution using the two convenience methods:
 756    using namespace boost::math;
 757    geometric g05d(0.5); // Using typedef - default type is double.
 758    geometric_distribution<> g05dd(0.5); // Using default RealType double.
 759
 760   // Basic sanity-check spot values.
 761
 762   // Test some simple double only examples.
 763   geometric_distribution<double> mydist(0.25);
 764   // success fraction == 0.25 == 25% or 1 in 4 successes.
 765   // Note: double values (matching the distribution definition) avoid the need for any casting.
 766
 767   // Check accessor functions return exact values for double at least.
 768   BOOST_CHECK_EQUAL(mydist.success_fraction(), static_cast<double>(1./4.));
 769
 770   //cout << numeric_limits<RealType>::epsilon() << endl;
 771
 772   // (Parameter value, arbitrarily zero, only communicates the floating point type).
 773 #ifdef TEST_FLOAT
 774   test_spots(0.0F); // Test float.
 775 #endif
 776 #ifdef TEST_DOUBLE
 777   test_spots(0.0); // Test double.
 778 #endif
 779 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 780 #ifdef TEST_LDOUBLE
 781   test_spots(0.0L); // Test long double.
 782 #endif
 783   #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
 784 #ifdef TEST_REAL_CONCEPT
 785     test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
 786 #endif
 787   #endif
 788 #else
 789    std::cout << "<note>The long double tests have been disabled on this platform "
 790       "either because the long double overloads of the usual math functions are "
 791       "not available at all, or because they are too inaccurate for these tests "
 792       "to pass.</note>" << std::endl;
 793 #endif
 794
 795
 796 } // BOOST_AUTO_TEST_CASE( test_main )
 797
 798 /*
 799
 800
 801
 802 */