libs/math/test/test_inverse_chi_squared_distribution.cpp

   1 // test_inverse_chi_squared.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 #ifdef _MSC_VER
  12 #  pragma warning (disable : 4310) // cast truncates constant value.
  13 #endif
  14
  15 // http://www.wolframalpha.com/input/?i=inverse+chisquare+distribution
  16
  17 #include <boost/math/tools/test.hpp>
  18 #include <boost/math/concepts/real_concept.hpp> // for real_concept
  19 using ::boost::math::concepts::real_concept;
  20
  21 //#include <boost/math/tools/test.hpp>
  22 #define BOOST_TEST_MAIN
  23 #include <boost/test/unit_test.hpp> // for test_main
  24 #include <boost/test/tools/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE_FRACTION
  25 #include "test_out_of_range.hpp"
  26
  27 #include <boost/math/distributions/inverse_chi_squared.hpp> // for inverse_chisquared_distribution
  28 using boost::math::inverse_chi_squared_distribution;
  29 using boost::math::cdf;
  30 using boost::math::pdf;
  31
  32 // Use Inverse Gamma distribution to check their relationship:
  33 // inverse_chi_squared<>(v) == inverse_gamma<>(v / 2., 0.5)
  34 #include <boost/math/distributions/inverse_gamma.hpp> // for inverse_gamma_distribution
  35 using boost::math::inverse_gamma_distribution;
  36 using boost::math::inverse_gamma;
  37 //  using  ::boost::math::cdf;
  38 //  using  ::boost::math::pdf;
  39
  40 #include <boost/math/special_functions/gamma.hpp>
  41 using boost::math::tgamma; // for naive pdf.
  42
  43 #include <iostream>
  44 using std::cout;
  45 using std::endl;
  46 #include <limits>
  47 using std::numeric_limits; // for epsilon.
  48
  49 template <class RealType>
  50 RealType naive_pdf(RealType df, RealType scale, RealType x)
  51 { // Formula from Wikipedia
  52    using namespace std; // For ADL of std functions.
  53    using boost::math::tgamma;
  54    RealType result = pow(scale * df/2, df/2) * exp(-df * scale/(2 * x));
  55    result /= tgamma(df/2) * pow(x, 1 + df/2);
  56    return result;
  57 }
  58
  59 // Test using a spot value from some other reference source,
  60 // in this case test values from output from R provided by Thomas Mang,
  61 // and Wolfram Mathematica by Mark Coleman.
  62
  63 template <class RealType>
  64 void test_spot(
  65      RealType degrees_of_freedom, // degrees_of_freedom,
  66      RealType scale, // scale,
  67      RealType x, // random variate x,
  68      RealType pd, // expected pdf,
  69      RealType P, // expected CDF,
  70      RealType Q, // expected complement of CDF,
  71      RealType tol) // test tolerance.
  72 {
  73    boost::math::inverse_chi_squared_distribution<RealType> dist(degrees_of_freedom, scale);
  74
  75    BOOST_CHECK_CLOSE_FRACTION
  76       ( // Compare to expected PDF.
  77       pdf(dist, x), // calculated.
  78       pd, // expected
  79       tol);
  80
  81    BOOST_CHECK_CLOSE_FRACTION( // Compare to naive pdf formula (probably less accurate).
  82       pdf(dist, x), naive_pdf(dist.degrees_of_freedom(), dist.scale(), x), tol);
  83
  84    BOOST_CHECK_CLOSE_FRACTION( // Compare to expected CDF.
  85       cdf(dist, x), P, tol);
  86
  87    if((P < 0.999) && (Q < 0.999))
  88    {  // We can only check this if P is not too close to 1,
  89       // so that we can guarantee Q is accurate:
  90       BOOST_CHECK_CLOSE_FRACTION(
  91         cdf(complement(dist, x)), Q, tol); // 1 - cdf
  92       BOOST_CHECK_CLOSE_FRACTION(
  93         quantile(dist, P), x, tol); // quantile(cdf) = x
  94       BOOST_CHECK_CLOSE_FRACTION(
  95         quantile(complement(dist, Q)), x, tol); // quantile(complement(1 - cdf)) = x
  96    }
  97 } // test_spot
  98
  99 template <class RealType> // Any floating-point type RealType.
 100 void test_spots(RealType)
 101 {
 102   // Basic sanity checks, some test data is to six decimal places only,
 103   // so set tolerance to 0.000001 (expressed as a percentage = 0.0001%).
 104
 105   RealType tolerance = 0.000001f;
 106   cout << "Tolerance = " << tolerance * 100 << "%." << endl;
 107
 108 // This test values from output from geoR (17 decimal digits) guided by Thomas Mang.
 109   test_spot(static_cast<RealType>(2), static_cast<RealType>(1./2.),
 110     // degrees_of_freedom, default scale = 1/df.
 111   static_cast<RealType>(1.L), // x.
 112   static_cast<RealType>(0.30326532985631671L), // pdf.
 113   static_cast<RealType>(0.60653065971263365L), // cdf.
 114   static_cast<RealType>(1 - 0.606530659712633657L), // cdf complement.
 115   tolerance  // tol
 116   );
 117
 118 // Tests from Mark Coleman & Georgi Boshnakov using Wolfram Mathematica.
 119   test_spot(static_cast<RealType>(10), static_cast<RealType>(0.1L), // degrees_of_freedom, scale
 120   static_cast<RealType>(0.2), // x
 121   static_cast<RealType>(1.6700235722635659824529759616528281217001163943570L), // pdf
 122   static_cast<RealType>(0.89117801891415124234834646836872197623907651175353L), // cdf
 123   static_cast<RealType>(1 - 0.89117801891415127L), // cdf complement
 124   tolerance  // tol
 125   );
 126
 127   test_spot(static_cast<RealType>(10), static_cast<RealType>(0.1L), // degrees_of_freedom, scale
 128   static_cast<RealType>(0.5), // x
 129   static_cast<RealType>(0.03065662009762021L), // pdf
 130   static_cast<RealType>(0.99634015317265628765454354418728984933240514654437L), // cdf
 131   static_cast<RealType>(1 - 0.99634015317265628765454354418728984933240514654437L), // cdf complement
 132   tolerance  // tol
 133   );
 134
 135
 136   test_spot(static_cast<RealType>(10), static_cast<RealType>(2), // degrees_of_freedom, scale
 137   static_cast<RealType>(0.5), // x
 138   static_cast<RealType>(0.00054964096598361569L), // pdf
 139   static_cast<RealType>(0.000016944743930067383903707995865261004246785511612700L), // cdf
 140   static_cast<RealType>(1 - 0.000016944743930067383903707995865261004246785511612700L), // cdf complement
 141   tolerance  // tol
 142   );
 143
 144   // Check some bad parameters to the distribution cause expected exception to be thrown.
 145 #ifndef BOOST_NO_EXCEPTIONS
 146   BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType> ichsqbad1(-1), std::domain_error); // negative degrees_of_freedom.
 147   BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType> ichsqbad2(1, -1), std::domain_error); // negative scale.
 148   BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType> ichsqbad3(-1, -1), std::domain_error); // negative scale and degrees_of_freedom.
 149 #else
 150   BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType>(-1), std::domain_error); // negative degrees_of_freedom.
 151   BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType>(1, -1), std::domain_error); // negative scale.
 152   BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType>(-1, -1), std::domain_error); // negative scale and degrees_of_freedom.
 153 #endif
 154   check_out_of_range<boost::math::inverse_chi_squared_distribution<RealType> >(1, 1);
 155
 156   inverse_chi_squared_distribution<RealType> ichsq;
 157
 158   if(std::numeric_limits<RealType>::has_infinity)
 159   {
 160     BOOST_MATH_CHECK_THROW(pdf(ichsq, +std::numeric_limits<RealType>::infinity()), std::domain_error); // x = + infinity, pdf = 0
 161     BOOST_MATH_CHECK_THROW(pdf(ichsq, -std::numeric_limits<RealType>::infinity()),  std::domain_error); // x = - infinity, pdf = 0
 162     BOOST_MATH_CHECK_THROW(cdf(ichsq, +std::numeric_limits<RealType>::infinity()),std::domain_error ); // x = + infinity, cdf = 1
 163     BOOST_MATH_CHECK_THROW(cdf(ichsq, -std::numeric_limits<RealType>::infinity()), std::domain_error); // x = - infinity, cdf = 0
 164     BOOST_MATH_CHECK_THROW(cdf(complement(ichsq, +std::numeric_limits<RealType>::infinity())), std::domain_error); // x = + infinity, c cdf = 0
 165     BOOST_MATH_CHECK_THROW(cdf(complement(ichsq, -std::numeric_limits<RealType>::infinity())), std::domain_error); // x = - infinity, c cdf = 1
 166 #ifndef BOOST_NO_EXCEPTIONS
 167     BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType> nbad1(std::numeric_limits<RealType>::infinity(), static_cast<RealType>(1)), std::domain_error); // +infinite mean
 168     BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType> nbad1(-std::numeric_limits<RealType>::infinity(),  static_cast<RealType>(1)), std::domain_error); // -infinite mean
 169     BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType> nbad1(static_cast<RealType>(0), std::numeric_limits<RealType>::infinity()), std::domain_error); // infinite sd
 170 #else
 171     BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType>(std::numeric_limits<RealType>::infinity(), static_cast<RealType>(1)), std::domain_error); // +infinite mean
 172     BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType>(-std::numeric_limits<RealType>::infinity(),  static_cast<RealType>(1)), std::domain_error); // -infinite mean
 173     BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType>(static_cast<RealType>(0), std::numeric_limits<RealType>::infinity()), std::domain_error); // infinite sd
 174 #endif
 175   }
 176
 177   if (std::numeric_limits<RealType>::has_quiet_NaN)
 178   { // If no longer allow x or p to be NaN, then these tests should throw.
 179     BOOST_MATH_CHECK_THROW(pdf(ichsq, +std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // x = NaN
 180     BOOST_MATH_CHECK_THROW(cdf(ichsq, +std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // x = NaN
 181     BOOST_MATH_CHECK_THROW(cdf(complement(ichsq, +std::numeric_limits<RealType>::quiet_NaN())), std::domain_error); // x = + infinity
 182     BOOST_MATH_CHECK_THROW(quantile(ichsq, std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // p = + quiet_NaN
 183     BOOST_MATH_CHECK_THROW(quantile(complement(ichsq, std::numeric_limits<RealType>::quiet_NaN())), std::domain_error); // p = + quiet_NaN
 184   }
 185     // Spot check for pdf using 'naive pdf' function
 186   for(RealType x = 0.5; x < 5; x += 0.5)
 187   {
 188     BOOST_CHECK_CLOSE_FRACTION(
 189       pdf(inverse_chi_squared_distribution<RealType>(5, 6), x),
 190       naive_pdf(RealType(5), RealType(6), x),
 191       tolerance);
 192   }   // Spot checks for parameters:
 193
 194   RealType tol_2eps = boost::math::tools::epsilon<RealType>() * 2; // 2 eps as a fraction.
 195   inverse_chi_squared_distribution<RealType> dist51(5, 1);
 196   inverse_chi_squared_distribution<RealType> dist52(5, 2);
 197   inverse_chi_squared_distribution<RealType> dist31(3, 1);
 198   inverse_chi_squared_distribution<RealType> dist111(11, 1);
 199   // 11 mean 0.10000000000000001, variance  0.0011111111111111111, sd 0.033333333333333333
 200
 201   using namespace std; // ADL of std names.
 202   using namespace boost::math;
 203
 204   inverse_chi_squared_distribution<RealType> dist10(10);
 205   //  mean, variance etc
 206   BOOST_CHECK_CLOSE_FRACTION(mean(dist10), static_cast<RealType>(0.125), tol_2eps);
 207   BOOST_CHECK_CLOSE_FRACTION(variance(dist10), static_cast<RealType>(0.0052083333333333333333333333333333333333333333333333L), tol_2eps);
 208   BOOST_CHECK_CLOSE_FRACTION(mode(dist10), static_cast<RealType>(0.08333333333333333333333333333333333333333333333L), tol_2eps);
 209   BOOST_CHECK_CLOSE_FRACTION(median(dist10), static_cast<RealType>(0.10704554778227709530244586234274024205738435512468L), tol_2eps);
 210   BOOST_CHECK_CLOSE_FRACTION(cdf(dist10, median(dist10)), static_cast<RealType>(0.5L), 4 * tol_2eps);
 211   BOOST_CHECK_CLOSE_FRACTION(skewness(dist10), static_cast<RealType>(3.4641016151377545870548926830117447338856105076208L), tol_2eps);
 212   BOOST_CHECK_CLOSE_FRACTION(kurtosis(dist10), static_cast<RealType>(45), tol_2eps);
 213   BOOST_CHECK_CLOSE_FRACTION(kurtosis_excess(dist10), static_cast<RealType>(45-3), tol_2eps);
 214
 215   tol_2eps = boost::math::tools::epsilon<RealType>() * 2; // 2 eps as a percentage.
 216
 217   // Special and limit cases:
 218
 219   RealType mx = (std::numeric_limits<RealType>::max)();
 220   RealType mi = (std::numeric_limits<RealType>::min)();
 221
 222   BOOST_CHECK_EQUAL(
 223   pdf(inverse_chi_squared_distribution<RealType>(1),
 224     static_cast<RealType>(mx)), // max()
 225     static_cast<RealType>(0)
 226     );
 227
 228   BOOST_CHECK_EQUAL(
 229   pdf(inverse_chi_squared_distribution<RealType>(1),
 230     static_cast<RealType>(mi)), // min()
 231     static_cast<RealType>(0)
 232     );
 233
 234   BOOST_CHECK_EQUAL(
 235     pdf(inverse_chi_squared_distribution<RealType>(1), static_cast<RealType>(0)), static_cast<RealType>(0));
 236   BOOST_CHECK_EQUAL(
 237     pdf(inverse_chi_squared_distribution<RealType>(3), static_cast<RealType>(0))
 238     , static_cast<RealType>(0.0f));
 239   BOOST_CHECK_EQUAL(
 240     cdf(inverse_chi_squared_distribution<RealType>(1), static_cast<RealType>(0))
 241     , static_cast<RealType>(0.0f));
 242   BOOST_CHECK_EQUAL(
 243     cdf(inverse_chi_squared_distribution<RealType>(2), static_cast<RealType>(0))
 244     , static_cast<RealType>(0.0f));
 245   BOOST_CHECK_EQUAL(
 246     cdf(inverse_chi_squared_distribution<RealType>(3L), static_cast<RealType>(0L))
 247     , static_cast<RealType>(0));
 248   BOOST_CHECK_EQUAL(
 249     cdf(complement(inverse_chi_squared_distribution<RealType>(1), static_cast<RealType>(0)))
 250     , static_cast<RealType>(1));
 251   BOOST_CHECK_EQUAL(
 252     cdf(complement(inverse_chi_squared_distribution<RealType>(2), static_cast<RealType>(0)))
 253     , static_cast<RealType>(1));
 254   BOOST_CHECK_EQUAL(
 255     cdf(complement(inverse_chi_squared_distribution<RealType>(3), static_cast<RealType>(0)))
 256     , static_cast<RealType>(1));
 257
 258   BOOST_MATH_CHECK_THROW(
 259     pdf(
 260     inverse_chi_squared_distribution<RealType>(static_cast<RealType>(-1)), // degrees_of_freedom negative.
 261     static_cast<RealType>(1)), std::domain_error
 262     );
 263   BOOST_MATH_CHECK_THROW(
 264     pdf(
 265     inverse_chi_squared_distribution<RealType>(static_cast<RealType>(8)),
 266     static_cast<RealType>(-1)), std::domain_error
 267     );
 268   BOOST_MATH_CHECK_THROW(
 269     cdf(
 270     inverse_chi_squared_distribution<RealType>(static_cast<RealType>(-1)),
 271     static_cast<RealType>(1)), std::domain_error
 272     );
 273   BOOST_MATH_CHECK_THROW(
 274     cdf(
 275     inverse_chi_squared_distribution<RealType>(static_cast<RealType>(8)),
 276     static_cast<RealType>(-1)), std::domain_error
 277     );
 278   BOOST_MATH_CHECK_THROW(
 279     cdf(complement(
 280     inverse_chi_squared_distribution<RealType>(static_cast<RealType>(-1)),
 281     static_cast<RealType>(1))), std::domain_error
 282     );
 283   BOOST_MATH_CHECK_THROW(
 284     cdf(complement(
 285     inverse_chi_squared_distribution<RealType>(static_cast<RealType>(8)),
 286     static_cast<RealType>(-1))), std::domain_error
 287     );
 288   BOOST_MATH_CHECK_THROW(
 289     quantile(
 290     inverse_chi_squared_distribution<RealType>(static_cast<RealType>(-1)),
 291     static_cast<RealType>(0.5)), std::domain_error
 292     );
 293   BOOST_MATH_CHECK_THROW(
 294     quantile(
 295     inverse_chi_squared_distribution<RealType>(static_cast<RealType>(8)),
 296     static_cast<RealType>(-1)), std::domain_error
 297     );
 298   BOOST_MATH_CHECK_THROW(
 299     quantile(
 300     inverse_chi_squared_distribution<RealType>(static_cast<RealType>(8)),
 301     static_cast<RealType>(1.1)), std::domain_error
 302     );
 303   BOOST_MATH_CHECK_THROW(
 304     quantile(complement(
 305     inverse_chi_squared_distribution<RealType>(static_cast<RealType>(-1)),
 306     static_cast<RealType>(0.5))), std::domain_error
 307     );
 308   BOOST_MATH_CHECK_THROW(
 309     quantile(complement(
 310     inverse_chi_squared_distribution<RealType>(static_cast<RealType>(8)),
 311     static_cast<RealType>(-1))), std::domain_error
 312     );
 313   BOOST_MATH_CHECK_THROW(
 314     quantile(complement(
 315     inverse_chi_squared_distribution<RealType>(static_cast<RealType>(8)),
 316     static_cast<RealType>(1.1))), std::domain_error
 317     );
 318 } // template <class RealType>void test_spots(RealType)
 319
 320
 321 BOOST_AUTO_TEST_CASE( test_main )
 322 {
 323   BOOST_MATH_CONTROL_FP;
 324
 325   double tol_few_eps = numeric_limits<double>::epsilon() * 4;
 326
 327   // Check that can generate inverse_chi_squared distribution using the two convenience methods:
 328   // inverse_chi_squared_distribution; // with default parameters, degrees_of_freedom = 1, scale - 1
 329   using boost::math::inverse_chi_squared;
 330
 331   // Some constructor tests using default double.
 332   double tol4eps = boost::math::tools::epsilon<double>() * 4; // 4 eps as a fraction.
 333
 334   inverse_chi_squared ichsqdef; // Using typedef and both default parameters.
 335
 336   BOOST_CHECK_EQUAL(ichsqdef.degrees_of_freedom(), 1.); // df == 1
 337   BOOST_CHECK_EQUAL(ichsqdef.scale(), 1); // scale == 1./df
 338   BOOST_CHECK_CLOSE_FRACTION(pdf(ichsqdef, 1), 0.24197072451914330, tol4eps);
 339   BOOST_CHECK_CLOSE_FRACTION(pdf(ichsqdef, 9), 0.013977156581221969, tol4eps);
 340
 341   inverse_chi_squared_distribution<double> ichisq102(10., 2); // Both parameters specified.
 342   BOOST_CHECK_EQUAL(ichisq102.degrees_of_freedom(), 10.); // Check both parameters stored OK.
 343   BOOST_CHECK_EQUAL(ichisq102.scale(), 2.); // Check both parameters stored OK.
 344
 345   inverse_chi_squared_distribution<double> ichisq10(10.); // Only df parameter specified (unscaled).
 346   BOOST_CHECK_EQUAL(ichisq10.degrees_of_freedom(), 10.); // Check  parameter stored.
 347   BOOST_CHECK_EQUAL(ichisq10.scale(), 0.1); // Check default scale = 1/df = 1/10 = 0.1
 348   BOOST_CHECK_CLOSE_FRACTION(pdf(ichisq10, 1),  0.00078975346316749169, tol4eps);
 349   BOOST_CHECK_CLOSE_FRACTION(pdf(ichisq10, 10), 0.0000000012385799798186384, tol4eps);
 350
 351   BOOST_CHECK_CLOSE_FRACTION(mode(ichisq10), 0.0833333333333333333333333333333333333333, tol4eps);
 352   // nu * xi / nu + 2 = 10 * 0.1 / (10 + 2) = 1/12 =  0.0833333...
 353   // mode is not defined in Mathematica.
 354   // See Discussion section http://en.wikipedia.org/wiki/Talk:Scaled-inverse-chi-square_distribution
 355   // for origin of this formula.
 356
 357   inverse_chi_squared_distribution<double> ichisq5(5.); // // Only df parameter specified.
 358   BOOST_CHECK_EQUAL(ichisq5.degrees_of_freedom(), 5.); // check  parameter stored.
 359   BOOST_CHECK_EQUAL(ichisq5.scale(), 1./5.); // check default is 1/df
 360   BOOST_CHECK_CLOSE_FRACTION(pdf(ichisq5, 0.2), 3.0510380337346841, tol4eps);
 361   BOOST_CHECK_CLOSE_FRACTION(cdf(ichisq5, 0.5), 0.84914503608460956, tol4eps);
 362   BOOST_CHECK_CLOSE_FRACTION(cdf(complement(ichisq5, 0.5)), 1 - 0.84914503608460956, tol4eps);
 363
 364   BOOST_CHECK_CLOSE_FRACTION(quantile(ichisq5, 0.84914503608460956), 0.5, tol4eps*100);
 365   BOOST_CHECK_CLOSE_FRACTION(quantile(complement(ichisq5, 1. - 0.84914503608460956)), 0.5, tol4eps*100);
 366
 367   // Check mean, etc spot values.
 368   inverse_chi_squared_distribution<double> ichisq81(8., 1.); // degrees_of_freedom = 5, scale = 1
 369   BOOST_CHECK_CLOSE_FRACTION(mean(ichisq81),1.33333333333333333333333333333333333333333, tol4eps);
 370   BOOST_CHECK_CLOSE_FRACTION(variance(ichisq81), 0.888888888888888888888888888888888888888888888, tol4eps);
 371   BOOST_CHECK_CLOSE_FRACTION(skewness(ichisq81), 2 * std::sqrt(8.), tol4eps);
 372   inverse_chi_squared_distribution<double> ichisq21(2., 1.);
 373   BOOST_CHECK_CLOSE_FRACTION(mode(ichisq21), 0.5, tol4eps);
 374   BOOST_CHECK_CLOSE_FRACTION(median(ichisq21), 1.4426950408889634, tol4eps);
 375
 376   inverse_chi_squared ichsq4(4.); // Using typedef and degrees_of_freedom parameter (and default scale = 1/df).
 377   BOOST_CHECK_EQUAL(ichsq4.degrees_of_freedom(), 4.); // df == 4.
 378   BOOST_CHECK_EQUAL(ichsq4.scale(), 0.25); // scale  == 1 /df == 1/4.
 379
 380   inverse_chi_squared ichsq32(3, 2);
 381   BOOST_CHECK_EQUAL(ichsq32.degrees_of_freedom(), 3.); // df == 3.
 382   BOOST_CHECK_EQUAL(ichsq32.scale(), 2); // scale  == 2
 383
 384   inverse_chi_squared ichsq11(1, 1); // Using explicit degrees_of_freedom parameter, and default scale = 1).
 385   BOOST_CHECK_CLOSE_FRACTION(mode(ichsq11), 0.3333333333333333333333333333333333333333, tol4eps);
 386   // (1 * 1)/ (1 + 2) = 1/3 using Wikipedia nu * xi /(nu + 2)
 387   BOOST_CHECK_EQUAL(ichsq11.degrees_of_freedom(), 1.); // df == 1 (default).
 388   BOOST_CHECK_EQUAL(ichsq11.scale(), 1.); // scale == 1.
 389   /*
 390   // Used to find some 'exact' values for testing mean, variance ...
 391   // First with scale fixed at unity (Wikipedia definition 1)
 392   cout << "df      scale            mean            variance              sd              median" << endl;
 393   for (int degrees_of_freedom = 8; degrees_of_freedom < 30; degrees_of_freedom++)
 394   {
 395     inverse_chi_squared ichisq(degrees_of_freedom, 1);
 396     cout.precision(17);
 397     cout << degrees_of_freedom << "    "  << 1 << "  " << mean(ichisq) << ' '
 398       << variance(ichisq) << ' ' << standard_deviation(ichisq)
 399       << ' ' << median(ichisq) << endl;
 400   }
 401
 402   // Default scale = 1 / df
 403   cout << "|\n" << "df           scale          mean            variance              sd              median" << endl;
 404   for (int degrees_of_freedom = 8; degrees_of_freedom < 30; degrees_of_freedom++)
 405   {
 406     inverse_chi_squared ichisq(degrees_of_freedom);
 407     cout.precision(17);
 408     cout << degrees_of_freedom << "    "  << 1./degrees_of_freedom << "  " << mean(ichisq) << ' '
 409       << variance(ichisq) << ' ' << standard_deviation(ichisq)
 410       << ' ' << median(ichisq) << endl;
 411   }
 412   */
 413   inverse_chi_squared_distribution<> ichisq14(14, 1); // Using default RealType double.
 414   BOOST_CHECK_CLOSE_FRACTION(mean(ichisq14), 1.166666666666666666666666666666666666666666666, tol4eps);
 415   BOOST_CHECK_CLOSE_FRACTION(variance(ichisq14), 0.272222222222222222222222222222222222222222222, tol4eps);
 416
 417   inverse_chi_squared_distribution<> ichisq121(12); // Using default RealType double.
 418   BOOST_CHECK_CLOSE_FRACTION(mean(ichisq121),  0.1, tol4eps);
 419   BOOST_CHECK_CLOSE_FRACTION(variance(ichisq121), 0.0025, tol4eps);
 420   BOOST_CHECK_CLOSE_FRACTION(standard_deviation(ichisq121), 0.05, tol4eps);
 421
 422   // and "using boost::math::inverse_chi_squared_distribution;".
 423   inverse_chi_squared_distribution<> ichsq23(2., 3.); // Using default RealType double.
 424   BOOST_CHECK_EQUAL(ichsq23.degrees_of_freedom(), 2.); //
 425   BOOST_CHECK_EQUAL(ichsq23.scale(), 3.); //
 426   BOOST_MATH_CHECK_THROW(mean(ichsq23), std::domain_error); // Degrees of freedom (nu) must be > 2
 427   BOOST_MATH_CHECK_THROW(variance(ichsq23), std::domain_error); // Degrees of freedom (nu) must be > 4
 428   BOOST_MATH_CHECK_THROW(skewness(ichsq23), std::domain_error); // Degrees of freedom (nu) must be > 6
 429   BOOST_MATH_CHECK_THROW(kurtosis_excess(ichsq23), std::domain_error); // Degrees of freedom (nu) must be > 8
 430
 431   { // Check relationship between inverse gamma and inverse chi_squared distributions.
 432   using boost::math::inverse_gamma_distribution;
 433
 434   double df = 2.;
 435   double scale = 1.;
 436   double alpha = df/2; // aka inv_gamma shape
 437   double beta = scale /2; // inv_gamma scale.
 438
 439   inverse_gamma_distribution<> ig(alpha, beta);
 440
 441   inverse_chi_squared_distribution<> ichsq(df, 1./df); // == default scale.
 442   BOOST_CHECK_EQUAL(pdf(ichsq, 0), 0); // Special case of zero x.
 443
 444   double x = 0.5;
 445   BOOST_CHECK_EQUAL(pdf(ig, x), pdf(ichsq, x)); // inv_gamma compared to inv_chisq
 446   BOOST_CHECK_EQUAL(cdf(ichsq, 0), 0); // Special case of zero.
 447   BOOST_CHECK_EQUAL(cdf(ig, x), cdf(ichsq, x)); // invgamma == invchisq
 448
 449   // Test pdf by comparing using naive_pdf with relation to inverse gamma distribution
 450   // wikipedia http://en.wikipedia.org/wiki/Scaled-inverse-chi-square_distribution related distributions.
 451   // So if naive_pdf is correct, inverse_chi_squared_distribution should agree.
 452   df = 1.; scale = 1.;
 453   BOOST_CHECK_CLOSE_FRACTION(naive_pdf(df, scale, x), pdf(ichsq11, x), tol_few_eps);
 454
 455   //inverse_gamma_distribution<> igd(df/2, (df * scale)/2);
 456   inverse_gamma_distribution<> igd11(df/2, df * scale/2);
 457   BOOST_CHECK_CLOSE_FRACTION(naive_pdf(df, scale, x), pdf(igd11, x), tol_few_eps);
 458   BOOST_CHECK_CLOSE_FRACTION(naive_pdf(df, scale, x), pdf(ichsq11, x), tol_few_eps);
 459
 460   df = 2; scale = 1;
 461   inverse_gamma_distribution<> igd21(df/2, df * scale/2);
 462   inverse_chi_squared_distribution<> ichsq21(df, scale);
 463   BOOST_CHECK_CLOSE_FRACTION(naive_pdf(df, scale, x), pdf(igd21, x), tol_few_eps); // 0.54134113294645081 OK
 464   BOOST_CHECK_CLOSE_FRACTION(naive_pdf(df, scale, x), pdf(ichsq21, x), tol_few_eps);
 465
 466   df = 2; scale = 2;
 467   inverse_gamma_distribution<> igd22(df/2, df * scale/2);
 468   inverse_chi_squared_distribution<> ichsq22(df, scale);
 469   BOOST_CHECK_CLOSE_FRACTION(naive_pdf(df, scale, x), pdf(igd22, x), tol_few_eps);
 470   BOOST_CHECK_CLOSE_FRACTION(naive_pdf(df, scale, x), pdf(ichsq22, x), tol_few_eps);
 471   }
 472
 473   // Check using float.
 474   inverse_chi_squared_distribution<float> igf23(1.f, 2.f); // Using explicit RealType float.
 475   BOOST_CHECK_EQUAL(igf23.degrees_of_freedom(), 1.f); //
 476   BOOST_CHECK_EQUAL(igf23.scale(), 2.f); //
 477
 478   // Check throws from bad parameters.
 479   inverse_chi_squared ig051(0.5, 1.); // degrees_of_freedom < 1, so wrong for mean.
 480   BOOST_MATH_CHECK_THROW(mean(ig051), std::domain_error);
 481   inverse_chi_squared ig191(1.9999, 1.); // degrees_of_freedom < 2, so wrong for variance.
 482   BOOST_MATH_CHECK_THROW(variance(ig191), std::domain_error);
 483   inverse_chi_squared ig291(2.9999, 1.); // degrees_of_freedom < 3, so wrong for skewness.
 484   BOOST_MATH_CHECK_THROW(skewness(ig291), std::domain_error);
 485   inverse_chi_squared ig391(3.9999, 1.); // degrees_of_freedom < 1, so wrong for kurtosis and kurtosis_excess.
 486   BOOST_MATH_CHECK_THROW(kurtosis(ig391), std::domain_error);
 487   BOOST_MATH_CHECK_THROW(kurtosis_excess(ig391), std::domain_error);
 488
 489   inverse_chi_squared ig102(10, 2); // Wolfram.com/ page 2, quantile = 2.96859.
 490   //http://reference.wolfram.com/mathematica/ref/InverseChiSquareDistribution.html
 491   BOOST_CHECK_CLOSE_FRACTION(quantile(ig102, 0.75), 2.96859, 0.000001);
 492   BOOST_CHECK_CLOSE_FRACTION(cdf(ig102, 2.96859), 0.75 , 0.000001);
 493   BOOST_CHECK_CLOSE_FRACTION(cdf(complement(ig102, 2.96859)), 1 - 0.75 , 0.00001);
 494   BOOST_CHECK_CLOSE_FRACTION(quantile(complement(ig102, 1 - 0.75)), 2.96859, 0.000001);
 495
 496   // Basic sanity-check spot values.
 497   // (Parameter value, arbitrarily zero, only communicates the floating point type).
 498   test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 %
 499   test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 %
 500 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
 501   test_spots(0.0L); // Test long double.
 502 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x0582))
 503   test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
 504 #endif
 505 #else
 506   std::cout << "<note>The long double tests have been disabled on this platform "
 507     "either because the long double overloads of the usual math functions are "
 508     "not available at all, or because they are too inaccurate for these tests "
 509     "to pass.</note>" << std::endl;
 510 #endif
 511
 512  /*    */
 513
 514 } // BOOST_AUTO_TEST_CASE( test_main )
 515
 516 /*
 517
 518 Output:
 519
 520
 521
 522
 523 */
 524
 525
 526