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

   1 // inverse_chi_squared_distribution_find_df_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 //#define BOOST_MATH_INSTRUMENT
  12
  13 // Example 1 of using inverse chi squared distribution
  14 #include <boost/math/distributions/inverse_chi_squared.hpp>
  15 using boost::math::inverse_chi_squared_distribution;  // inverse_chi_squared_distribution.
  16 using boost::math::inverse_chi_squared; //typedef for nverse_chi_squared_distribution double.
  17
  18 #include <iostream>
  19 using std::cout;    using std::endl;
  20 #include <iomanip>
  21 using std::setprecision;
  22 using std::setw;
  23 #include <cmath>
  24 using std::sqrt;
  25
  26 int main()
  27 {
  28   cout << "Example using Inverse chi squared distribution to find df. " << endl;
  29   try
  30   {
  31     cout.precision(std::numeric_limits<double>::max_digits10); //
  32     int i = std::numeric_limits<double>::max_digits10;
  33     cout << "Show all potentially significant decimal digits std::numeric_limits<double>::max_digits10 = " << i << endl;
  34
  35     cout.precision(3);
  36     double nu = 10.;
  37     double scale1 = 1./ nu; // 1st definition sigma^2 = 1/df;
  38     double scale2 = 1.; // 2nd definition sigma^2 = 1
  39     inverse_chi_squared sichsq(nu, 1/nu); // Explicitly scaled to default scale = 1/df.
  40     inverse_chi_squared ichsq(nu); // Implicitly scaled to default scale = 1/df.
  41     // Try degrees of freedom estimator
  42
  43     //double df = chi_squared::find_degrees_of_freedom(-diff, alpha[i], alpha[i], variance);
  44
  45     cout << "ichsq.degrees_of_freedom() = " << ichsq.degrees_of_freedom() << endl;
  46
  47     double diff = 0.5;  // difference from variance to detect (delta).
  48     double variance = 1.; // true variance
  49     double alpha = 0.9;
  50     double beta = 0.9;
  51
  52     cout << "diff = " << diff
  53       << ", variance = " << variance << ", ratio = " << diff/variance
  54       << ", alpha = " << alpha << ", beta = " << beta << endl;
  55
  56     /* inverse_chi_square_df_estimator is not in the code base anymore?
  57
  58     using boost::math::detail::inverse_chi_square_df_estimator;
  59     using boost::math::policies::default_policy;
  60     inverse_chi_square_df_estimator<> a_df(alpha, beta, variance, diff);
  61
  62     cout << "df    est" << endl;
  63     for (double df = 1; df < 3; df += 0.1)
  64     {
  65       double est_df = a_df(1);
  66       cout << df << "    " << a_df(df) << endl;
  67     }
  68     */
  69
  70     //template <class F, class T, class Tol, class Policy>std::pair<T, T>
  71     // bracket_and_solve_root(F f, const T& guess, T factor, bool rising, Tol tol, std::uintmax_t& max_iter, const Policy& pol)
  72
  73
  74     // TODO: Not implemented
  75     //double df = inverse_chi_squared_distribution<>::find_degrees_of_freedom(diff, alpha, beta, variance, 0);
  76     //cout << df << endl;
  77   }
  78   catch(const std::exception& e)
  79   { // Always useful to include try & catch blocks because default policies
  80     // are to throw exceptions on arguments that cause errors like underflow, overflow.
  81     // Lacking try & catch blocks, the program will abort without a message below,
  82     // which may give some helpful clues as to the cause of the exception.
  83     std::cout <<
  84       "\n""Message from thrown exception was:\n   " << e.what() << std::endl;
  85   }
  86   return 0;
  87 }  // int main()
  88
  89 /*
  90
  91 Output is:
  92
  93   Example using Inverse chi squared distribution to find df.
  94   Show all potentially significant decimal digits std::numeric_limits<double>::max_digits10 = 17
  95   10
  96
  97   Message from thrown exception was:
  98      Error in function boost::math::inverse_chi_squared_distribution<double>::inverse_chi_squared_distribution: Degrees of freedom argument is 1.#INF, but must be > 0 !
  99 diff = 0.5, variance = 1, ratio = 0.5, alpha = 0.1, beta = 0.1
 100   df    est
 101   1    1
 102   ratio+1 = 1.5, quantile(0.1) = 0.00618, cdf = 6.5e-037, result = -0.1
 103   1.1    -0.1
 104   ratio+1 = 1.5, quantile(0.1) = 0.00903, cdf = 1.2e-025, result = -0.1
 105   1.2    -0.1
 106   ratio+1 = 1.5, quantile(0.1) = 0.0125, cdf = 8.25e-019, result = -0.1
 107   1.3    -0.1
 108   ratio+1 = 1.5, quantile(0.1) = 0.0166, cdf = 2.17e-014, result = -0.1
 109   1.4    -0.1
 110   ratio+1 = 1.5, quantile(0.1) = 0.0212, cdf = 2.2e-011, result = -0.1
 111   1.5    -0.1
 112   ratio+1 = 1.5, quantile(0.1) = 0.0265, cdf = 3e-009, result = -0.1
 113   1.6    -0.1
 114   ratio+1 = 1.5, quantile(0.1) = 0.0323, cdf = 1.11e-007, result = -0.1
 115   1.7    -0.1
 116   ratio+1 = 1.5, quantile(0.1) = 0.0386, cdf = 1.7e-006, result = -0.1
 117   1.8    -0.1
 118   ratio+1 = 1.5, quantile(0.1) = 0.0454, cdf = 1.41e-005, result = -0.1
 119   1.9    -0.1
 120   ratio+1 = 1.5, quantile(0.1) = 0.0527, cdf = 7.55e-005, result = -0.1
 121   2    -0.1
 122   ratio+1 = 1.5, quantile(0.1) = 0.0604, cdf = 0.000291, result = -0.1
 123   2.1    -0.1
 124   ratio+1 = 1.5, quantile(0.1) = 0.0685, cdf = 0.00088, result = -0.1
 125   2.2    -0.1
 126   ratio+1 = 1.5, quantile(0.1) = 0.0771, cdf = 0.0022, result = -0.0999
 127   2.3    -0.0999
 128   ratio+1 = 1.5, quantile(0.1) = 0.0859, cdf = 0.00475, result = -0.0997
 129   2.4    -0.0997
 130   ratio+1 = 1.5, quantile(0.1) = 0.0952, cdf = 0.00911, result = -0.0993
 131   2.5    -0.0993
 132   ratio+1 = 1.5, quantile(0.1) = 0.105, cdf = 0.0159, result = -0.0984
 133   2.6    -0.0984
 134   ratio+1 = 1.5, quantile(0.1) = 0.115, cdf = 0.0257, result = -0.0967
 135   2.7    -0.0967
 136   ratio+1 = 1.5, quantile(0.1) = 0.125, cdf = 0.039, result = -0.094
 137   2.8    -0.094
 138   ratio+1 = 1.5, quantile(0.1) = 0.135, cdf = 0.056, result = -0.0897
 139   2.9    -0.0897
 140   ratio+1 = 1.5, quantile(0.1) = 20.6, cdf = 1, result = 0.9
 141
 142     ichsq.degrees_of_freedom() = 10
 143   diff = 0.5, variance = 1, ratio = 0.5, alpha = 0.9, beta = 0.9
 144   df    est
 145   1    1
 146   ratio+1 = 1.5, quantile(0.9) = 0.729, cdf = 0.269, result = -0.729
 147   1.1    -0.729
 148   ratio+1 = 1.5, quantile(0.9) = 0.78, cdf = 0.314, result = -0.693
 149   1.2    -0.693
 150   ratio+1 = 1.5, quantile(0.9) = 0.83, cdf = 0.36, result = -0.655
 151   1.3    -0.655
 152   ratio+1 = 1.5, quantile(0.9) = 0.879, cdf = 0.405, result = -0.615
 153   1.4    -0.615
 154   ratio+1 = 1.5, quantile(0.9) = 0.926, cdf = 0.449, result = -0.575
 155   1.5    -0.575
 156   ratio+1 = 1.5, quantile(0.9) = 0.973, cdf = 0.492, result = -0.535
 157   1.6    -0.535
 158   ratio+1 = 1.5, quantile(0.9) = 1.02, cdf = 0.534, result = -0.495
 159   1.7    -0.495
 160   ratio+1 = 1.5, quantile(0.9) = 1.06, cdf = 0.574, result = -0.455
 161   1.8    -0.455
 162   ratio+1 = 1.5, quantile(0.9) = 1.11, cdf = 0.612, result = -0.417
 163   1.9    -0.417
 164   ratio+1 = 1.5, quantile(0.9) = 1.15, cdf = 0.648, result = -0.379
 165   2    -0.379
 166   ratio+1 = 1.5, quantile(0.9) = 1.19, cdf = 0.681, result = -0.342
 167   2.1    -0.342
 168   ratio+1 = 1.5, quantile(0.9) = 1.24, cdf = 0.713, result = -0.307
 169   2.2    -0.307
 170   ratio+1 = 1.5, quantile(0.9) = 1.28, cdf = 0.742, result = -0.274
 171   2.3    -0.274
 172   ratio+1 = 1.5, quantile(0.9) = 1.32, cdf = 0.769, result = -0.242
 173   2.4    -0.242
 174   ratio+1 = 1.5, quantile(0.9) = 1.36, cdf = 0.793, result = -0.212
 175   2.5    -0.212
 176   ratio+1 = 1.5, quantile(0.9) = 1.4, cdf = 0.816, result = -0.184
 177   2.6    -0.184
 178   ratio+1 = 1.5, quantile(0.9) = 1.44, cdf = 0.836, result = -0.157
 179   2.7    -0.157
 180   ratio+1 = 1.5, quantile(0.9) = 1.48, cdf = 0.855, result = -0.133
 181   2.8    -0.133
 182   ratio+1 = 1.5, quantile(0.9) = 1.52, cdf = 0.872, result = -0.11
 183   2.9    -0.11
 184   ratio+1 = 1.5, quantile(0.9) = 29.6, cdf = 1, result = 0.1
 185
 186
 187   */
 188