ceph/src/boost/libs/math/example/inverse_chi_squared_example.cpp

   1 // inverse_chi_squared_distribution_example.cpp
   2
   3 // Copyright Paul A. Bristow 2010.
   4 // Copyright Thomas Mang 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 // Example 1 of using inverse chi squared distribution
  12 #include <boost/math/distributions/inverse_chi_squared.hpp>
  13 using boost::math::inverse_chi_squared_distribution;  // inverse_chi_squared_distribution.
  14 using boost::math::inverse_chi_squared; //typedef for nverse_chi_squared_distribution double.
  15
  16 #include <iostream>
  17 using std::cout;    using std::endl;
  18 #include <iomanip>
  19 using std::setprecision;
  20 using std::setw;
  21 #include <cmath>
  22 using std::sqrt;
  23
  24 template <class RealType>
  25 RealType naive_pdf1(RealType df, RealType x)
  26 { // Formula from Wikipedia http://en.wikipedia.org/wiki/Inverse-chi-square_distribution
  27   // definition 1 using tgamma for simplicity as a check.
  28    using namespace std; // For ADL of std functions.
  29    using boost::math::tgamma;
  30    RealType df2 = df / 2;
  31    RealType result = (pow(2., -df2) * pow(x, (-df2 -1)) * exp(-1/(2 * x) ) )
  32       / tgamma(df2);  //
  33    return result;
  34 }
  35
  36 template <class RealType>
  37 RealType naive_pdf2(RealType df, RealType x)
  38 { // Formula from Wikipedia http://en.wikipedia.org/wiki/Inverse-chi-square_distribution
  39   // Definition 2, using tgamma for simplicity as a check.
  40   // Not scaled, so assumes scale = 1 and only uses nu aka df.
  41    using namespace std; // For ADL of std functions.
  42    using boost::math::tgamma;
  43    RealType df2 = df / 2;
  44    RealType result = (pow(df2, df2) * pow(x, (-df2 -1)) * exp(-df/(2*x) ) )
  45      / tgamma(df2);
  46    return result;
  47 }
  48
  49 template <class RealType>
  50 RealType naive_pdf3(RealType df, RealType scale, RealType x)
  51 { // Formula from Wikipedia http://en.wikipedia.org/wiki/Scaled-inverse-chi-square_distribution
  52   // *Scaled* version, definition 3, df aka nu, scale aka sigma^2
  53   // using tgamma for simplicity as a check.
  54    using namespace std; // For ADL of std functions.
  55    using boost::math::tgamma;
  56    RealType df2 = df / 2;
  57    RealType result = (pow(scale * df2, df2) * exp(-df2 * scale/x) )
  58      / (tgamma(df2) * pow(x, 1+df2));
  59    return result;
  60 }
  61
  62 template <class RealType>
  63 RealType naive_pdf4(RealType df, RealType scale, RealType x)
  64 { // Formula from http://mathworld.wolfram.com/InverseChi-SquaredDistribution.html
  65   // Weisstein, Eric W. "Inverse Chi-Squared Distribution." From MathWorld--A Wolfram Web Resource.
  66   // *Scaled* version, definition 3, df aka nu, scale aka sigma^2
  67   // using tgamma for simplicity as a check.
  68    using namespace std; // For ADL of std functions.
  69    using boost::math::tgamma;
  70    RealType nu = df; // Wolfram uses greek symbols nu,
  71    RealType xi = scale; // and xi.
  72    RealType result =
  73      pow(2, -nu/2) *  exp(- (nu * xi)/(2 * x)) * pow(x, -1-nu/2) * pow((nu * xi), nu/2)
  74      / tgamma(nu/2);
  75    return result;
  76 }
  77
  78 int main()
  79 {
  80
  81   cout << "Example (basic) using Inverse chi squared distribution. " << endl;
  82
  83   // TODO find a more practical/useful example.  Suggestions welcome?
  84
  85 #ifdef BOOST_NO_CXX11_NUMERIC_LIMITS
  86   int max_digits10 = 2 + (boost::math::policies::digits<double, boost::math::policies::policy<> >() * 30103UL) / 100000UL;
  87   cout << "BOOST_NO_CXX11_NUMERIC_LIMITS is defined" << endl;
  88 #else
  89   int max_digits10 = std::numeric_limits<double>::max_digits10;
  90 #endif
  91   cout << "Show all potentially significant decimal digits std::numeric_limits<double>::max_digits10 = "
  92     << max_digits10 << endl;
  93   cout.precision(max_digits10); //
  94
  95   inverse_chi_squared ichsqdef; // All defaults  - not very useful!
  96   cout << "default df = " << ichsqdef.degrees_of_freedom()
  97     << ", default scale = " <<  ichsqdef.scale() << endl; //  default df = 1, default scale = 0.5
  98
  99    inverse_chi_squared ichsqdef4(4); // Unscaled version, default scale = 1 / degrees_of_freedom
 100    cout << "default df = " << ichsqdef4.degrees_of_freedom()
 101     << ", default scale = " <<  ichsqdef4.scale() << endl; //  default df = 4, default scale = 2
 102
 103    inverse_chi_squared ichsqdef32(3, 2); // Scaled version, both degrees_of_freedom and scale specified.
 104    cout << "default df = " << ichsqdef32.degrees_of_freedom()
 105     << ", default scale = " <<  ichsqdef32.scale() << endl; // default df = 3, default scale = 2
 106
 107   {
 108     cout.precision(3);
 109     double nu = 5.;
 110     //double scale1 = 1./ nu; // 1st definition sigma^2 = 1/df;
 111     //double scale2 = 1.; // 2nd definition sigma^2 = 1
 112     inverse_chi_squared ichsq(nu, 1/nu); // Not scaled
 113     inverse_chi_squared sichsq(nu, 1/nu); // scaled
 114
 115     cout << "nu = " << ichsq.degrees_of_freedom() << ", scale = " << ichsq.scale() << endl;
 116
 117     int width = 8;
 118
 119     cout << "     x        pdf      pdf1   pdf2  pdf(scaled)    pdf       pdf      cdf     cdf" << endl;
 120     for (double x = 0.0; x < 1.; x += 0.1)
 121     {
 122       cout
 123         << setw(width) << x
 124         << ' ' << setw(width) << pdf(ichsq, x) // unscaled
 125         << ' ' << setw(width) << naive_pdf1(nu,  x) // Wiki def 1 unscaled matches graph
 126         << ' ' << setw(width) << naive_pdf2(nu,  x) // scale = 1 - 2nd definition.
 127         << ' ' << setw(width) << naive_pdf3(nu, 1/nu, x) // scaled
 128         << ' ' << setw(width) << naive_pdf4(nu, 1/nu, x) // scaled
 129         << ' ' << setw(width) << pdf(sichsq, x)  // scaled
 130         << ' ' << setw(width) << cdf(sichsq, x)  // scaled
 131         << ' ' << setw(width) << cdf(ichsq, x)  // unscaled
 132        << endl;
 133     }
 134   }
 135
 136   cout.precision(max_digits10);
 137
 138   inverse_chi_squared ichisq(2., 0.5);
 139   cout << "pdf(ichisq, 1.) = " << pdf(ichisq, 1.) << endl;
 140   cout << "cdf(ichisq, 1.) = " << cdf(ichisq, 1.) << endl;
 141
 142   return 0;
 143 }  // int main()
 144
 145 /*
 146
 147 Output is:
 148  Example (basic) using Inverse chi squared distribution.
 149   Show all potentially significant decimal digits std::numeric_limits<double>::max_digits10 = 17
 150   default df = 1, default scale = 1
 151   default df = 4, default scale = 0.25
 152   default df = 3, default scale = 2
 153   nu = 5, scale = 0.2
 154        x        pdf      pdf1   pdf2  pdf(scaled)    pdf       pdf      cdf     cdf
 155          0        0    -1.#J    -1.#J    -1.#J    -1.#J        0        0        0
 156        0.1     2.83     2.83 3.26e-007     2.83     2.83     2.83   0.0752   0.0752
 157        0.2     3.05     3.05  0.00774     3.05     3.05     3.05    0.416    0.416
 158        0.3      1.7      1.7    0.121      1.7      1.7      1.7    0.649    0.649
 159        0.4    0.941    0.941    0.355    0.941    0.941    0.941    0.776    0.776
 160        0.5    0.553    0.553    0.567    0.553    0.553    0.553    0.849    0.849
 161        0.6    0.345    0.345    0.689    0.345    0.345    0.345    0.893    0.893
 162        0.7    0.227    0.227    0.728    0.227    0.227    0.227    0.921    0.921
 163        0.8    0.155    0.155    0.713    0.155    0.155    0.155     0.94     0.94
 164        0.9     0.11     0.11    0.668     0.11     0.11     0.11    0.953    0.953
 165          1   0.0807   0.0807     0.61   0.0807   0.0807   0.0807    0.963    0.963
 166   pdf(ichisq, 1.) = 0.30326532985631671
 167   cdf(ichisq, 1.) = 0.60653065971263365
 168
 169
 170 */
 171