ceph/src/boost/boost/math/distributions/inverse_chi_squared.hpp

   1 // Copyright John Maddock 2010.
   2 // Copyright Paul A. Bristow 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 #ifndef BOOST_MATH_DISTRIBUTIONS_INVERSE_CHI_SQUARED_HPP
  10 #define BOOST_MATH_DISTRIBUTIONS_INVERSE_CHI_SQUARED_HPP
  11
  12 #include <boost/math/distributions/fwd.hpp>
  13 #include <boost/math/special_functions/gamma.hpp> // for incomplete beta.
  14 #include <boost/math/distributions/complement.hpp> // for complements.
  15 #include <boost/math/distributions/detail/common_error_handling.hpp> // for error checks.
  16 #include <boost/math/special_functions/fpclassify.hpp> // for isfinite
  17
  18 // See http://en.wikipedia.org/wiki/Scaled-inverse-chi-square_distribution
  19 // for definitions of this scaled version.
  20 // See http://en.wikipedia.org/wiki/Inverse-chi-square_distribution
  21 // for unscaled version.
  22
  23 // http://reference.wolfram.com/mathematica/ref/InverseChiSquareDistribution.html
  24 // Weisstein, Eric W. "Inverse Chi-Squared Distribution." From MathWorld--A Wolfram Web Resource.
  25 // http://mathworld.wolfram.com/InverseChi-SquaredDistribution.html
  26
  27 #include <utility>
  28
  29 namespace boost{ namespace math{
  30
  31 namespace detail
  32 {
  33   template <class RealType, class Policy>
  34   inline bool check_inverse_chi_squared( // Check both distribution parameters.
  35         const char* function,
  36         RealType degrees_of_freedom, // degrees_of_freedom (aka nu).
  37         RealType scale,  // scale (aka sigma^2)
  38         RealType* result,
  39         const Policy& pol)
  40   {
  41      return check_scale(function, scale, result, pol)
  42        && check_df(function, degrees_of_freedom,
  43        result, pol);
  44   } // bool check_inverse_chi_squared
  45 } // namespace detail
  46
  47 template <class RealType = double, class Policy = policies::policy<> >
  48 class inverse_chi_squared_distribution
  49 {
  50 public:
  51    typedef RealType value_type;
  52    typedef Policy policy_type;
  53
  54    inverse_chi_squared_distribution(RealType df, RealType l_scale) : m_df(df), m_scale (l_scale)
  55    {
  56       RealType result;
  57       detail::check_df(
  58          "boost::math::inverse_chi_squared_distribution<%1%>::inverse_chi_squared_distribution",
  59          m_df, &result, Policy())
  60          && detail::check_scale(
  61 "boost::math::inverse_chi_squared_distribution<%1%>::inverse_chi_squared_distribution",
  62          m_scale, &result,  Policy());
  63    } // inverse_chi_squared_distribution constructor
  64
  65    inverse_chi_squared_distribution(RealType df = 1) : m_df(df)
  66    {
  67       RealType result;
  68       m_scale = 1 / m_df ; // Default scale = 1 / degrees of freedom (Wikipedia definition 1).
  69       detail::check_df(
  70          "boost::math::inverse_chi_squared_distribution<%1%>::inverse_chi_squared_distribution",
  71          m_df, &result, Policy());
  72    } // inverse_chi_squared_distribution
  73
  74    RealType degrees_of_freedom()const
  75    {
  76       return m_df; // aka nu
  77    }
  78    RealType scale()const
  79    {
  80       return m_scale;  // aka xi
  81    }
  82
  83    // Parameter estimation:  NOT implemented yet.
  84    //static RealType find_degrees_of_freedom(
  85    //   RealType difference_from_variance,
  86    //   RealType alpha,
  87    //   RealType beta,
  88    //   RealType variance,
  89    //   RealType hint = 100);
  90
  91 private:
  92    // Data members:
  93    RealType m_df;  // degrees of freedom are treated as a real number.
  94    RealType m_scale;  // distribution scale.
  95
  96 }; // class chi_squared_distribution
  97
  98 typedef inverse_chi_squared_distribution<double> inverse_chi_squared;
  99
 100 #ifdef __cpp_deduction_guides
 101 template <class RealType>
 102 inverse_chi_squared_distribution(RealType)->inverse_chi_squared_distribution<typename boost::math::tools::promote_args<RealType>::type>;
 103 template <class RealType>
 104 inverse_chi_squared_distribution(RealType,RealType)->inverse_chi_squared_distribution<typename boost::math::tools::promote_args<RealType>::type>;
 105 #endif
 106
 107 template <class RealType, class Policy>
 108 inline const std::pair<RealType, RealType> range(const inverse_chi_squared_distribution<RealType, Policy>& /*dist*/)
 109 {  // Range of permissible values for random variable x.
 110    using boost::math::tools::max_value;
 111    return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // 0 to + infinity.
 112 }
 113
 114 template <class RealType, class Policy>
 115 inline const std::pair<RealType, RealType> support(const inverse_chi_squared_distribution<RealType, Policy>& /*dist*/)
 116 {  // Range of supported values for random variable x.
 117    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
 118    return std::pair<RealType, RealType>(static_cast<RealType>(0), tools::max_value<RealType>()); // 0 to + infinity.
 119 }
 120
 121 template <class RealType, class Policy>
 122 RealType pdf(const inverse_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
 123 {
 124    BOOST_MATH_STD_USING  // for ADL of std functions.
 125    RealType df = dist.degrees_of_freedom();
 126    RealType scale = dist.scale();
 127    RealType error_result;
 128
 129    static const char* function = "boost::math::pdf(const inverse_chi_squared_distribution<%1%>&, %1%)";
 130
 131    if(false == detail::check_inverse_chi_squared
 132      (function, df, scale, &error_result, Policy())
 133      )
 134    { // Bad distribution.
 135       return error_result;
 136    }
 137    if((x < 0) || !(boost::math::isfinite)(x))
 138    { // Bad x.
 139       return policies::raise_domain_error<RealType>(
 140          function, "inverse Chi Square parameter was %1%, but must be >= 0 !", x, Policy());
 141    }
 142
 143    if(x == 0)
 144    { // Treat as special case.
 145      return 0;
 146    }
 147    // Wikipedia scaled inverse chi sq (df, scale) related to inv gamma (df/2, df * scale /2)
 148    // so use inverse gamma pdf with shape = df/2, scale df * scale /2
 149    // RealType shape = df /2; // inv_gamma shape
 150    // RealType scale = df * scale/2; // inv_gamma scale
 151    // RealType result = gamma_p_derivative(shape, scale / x, Policy()) * scale / (x * x);
 152    RealType result = df * scale/2 / x;
 153    if(result < tools::min_value<RealType>())
 154       return 0; // Random variable is near enough infinite.
 155    result = gamma_p_derivative(df/2, result, Policy()) * df * scale/2;
 156    if(result != 0) // prevent 0 / 0,  gamma_p_derivative -> 0 faster than x^2
 157       result /= (x * x);
 158    return result;
 159 } // pdf
 160
 161 template <class RealType, class Policy>
 162 inline RealType cdf(const inverse_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
 163 {
 164    static const char* function = "boost::math::cdf(const inverse_chi_squared_distribution<%1%>&, %1%)";
 165    RealType df = dist.degrees_of_freedom();
 166    RealType scale = dist.scale();
 167    RealType error_result;
 168
 169    if(false ==
 170        detail::check_inverse_chi_squared(function, df, scale, &error_result, Policy())
 171      )
 172    { // Bad distribution.
 173       return error_result;
 174    }
 175    if((x < 0) || !(boost::math::isfinite)(x))
 176    { // Bad x.
 177       return policies::raise_domain_error<RealType>(
 178          function, "inverse Chi Square parameter was %1%, but must be >= 0 !", x, Policy());
 179    }
 180    if (x == 0)
 181    { // Treat zero as a special case.
 182      return 0;
 183    }
 184    // RealType shape = df /2; // inv_gamma shape,
 185    // RealType scale = df * scale/2; // inv_gamma scale,
 186    // result = boost::math::gamma_q(shape, scale / x, Policy()); // inverse_gamma code.
 187    return boost::math::gamma_q(df / 2, (df * (scale / 2)) / x, Policy());
 188 } // cdf
 189
 190 template <class RealType, class Policy>
 191 inline RealType quantile(const inverse_chi_squared_distribution<RealType, Policy>& dist, const RealType& p)
 192 {
 193    using boost::math::gamma_q_inv;
 194    RealType df = dist.degrees_of_freedom();
 195    RealType scale = dist.scale();
 196
 197    static const char* function = "boost::math::quantile(const inverse_chi_squared_distribution<%1%>&, %1%)";
 198    // Error check:
 199    RealType error_result;
 200    if(false == detail::check_df(
 201          function, df, &error_result, Policy())
 202          && detail::check_probability(
 203             function, p, &error_result, Policy()))
 204    {
 205       return error_result;
 206    }
 207    if(false == detail::check_probability(
 208             function, p, &error_result, Policy()))
 209    {
 210       return error_result;
 211    }
 212    // RealType shape = df /2; // inv_gamma shape,
 213    // RealType scale = df * scale/2; // inv_gamma scale,
 214    // result = scale / gamma_q_inv(shape, p, Policy());
 215       RealType result = gamma_q_inv(df /2, p, Policy());
 216       if(result == 0)
 217          return policies::raise_overflow_error<RealType, Policy>(function, "Random variable is infinite.", Policy());
 218       result = df * (scale / 2) / result;
 219       return result;
 220 } // quantile
 221
 222 template <class RealType, class Policy>
 223 inline RealType cdf(const complemented2_type<inverse_chi_squared_distribution<RealType, Policy>, RealType>& c)
 224 {
 225    using boost::math::gamma_q_inv;
 226    RealType const& df = c.dist.degrees_of_freedom();
 227    RealType const& scale = c.dist.scale();
 228    RealType const& x = c.param;
 229    static const char* function = "boost::math::cdf(const inverse_chi_squared_distribution<%1%>&, %1%)";
 230    // Error check:
 231    RealType error_result;
 232    if(false == detail::check_df(
 233          function, df, &error_result, Policy()))
 234    {
 235       return error_result;
 236    }
 237    if (x == 0)
 238    { // Treat zero as a special case.
 239      return 1;
 240    }
 241    if((x < 0) || !(boost::math::isfinite)(x))
 242    {
 243       return policies::raise_domain_error<RealType>(
 244          function, "inverse Chi Square parameter was %1%, but must be > 0 !", x, Policy());
 245    }
 246    // RealType shape = df /2; // inv_gamma shape,
 247    // RealType scale = df * scale/2; // inv_gamma scale,
 248    // result = gamma_p(shape, scale/c.param, Policy()); use inv_gamma.
 249
 250    return gamma_p(df / 2, (df * scale/2) / x, Policy()); // OK
 251 } // cdf(complemented
 252
 253 template <class RealType, class Policy>
 254 inline RealType quantile(const complemented2_type<inverse_chi_squared_distribution<RealType, Policy>, RealType>& c)
 255 {
 256    using boost::math::gamma_q_inv;
 257
 258    RealType const& df = c.dist.degrees_of_freedom();
 259    RealType const& scale = c.dist.scale();
 260    RealType const& q = c.param;
 261    static const char* function = "boost::math::quantile(const inverse_chi_squared_distribution<%1%>&, %1%)";
 262    // Error check:
 263    RealType error_result;
 264    if(false == detail::check_df(function, df, &error_result, Policy()))
 265    {
 266       return error_result;
 267    }
 268    if(false == detail::check_probability(function, q, &error_result, Policy()))
 269    {
 270       return error_result;
 271    }
 272    // RealType shape = df /2; // inv_gamma shape,
 273    // RealType scale = df * scale/2; // inv_gamma scale,
 274    // result = scale / gamma_p_inv(shape, q, Policy());  // using inv_gamma.
 275    RealType result = gamma_p_inv(df/2, q, Policy());
 276    if(result == 0)
 277       return policies::raise_overflow_error<RealType, Policy>(function, "Random variable is infinite.", Policy());
 278    result = (df * scale / 2) / result;
 279    return result;
 280 } // quantile(const complement
 281
 282 template <class RealType, class Policy>
 283 inline RealType mean(const inverse_chi_squared_distribution<RealType, Policy>& dist)
 284 { // Mean of inverse Chi-Squared distribution.
 285    RealType df = dist.degrees_of_freedom();
 286    RealType scale = dist.scale();
 287
 288    static const char* function = "boost::math::mean(const inverse_chi_squared_distribution<%1%>&)";
 289    if(df <= 2)
 290       return policies::raise_domain_error<RealType>(
 291          function,
 292          "inverse Chi-Squared distribution only has a mode for degrees of freedom > 2, but got degrees of freedom = %1%.",
 293          df, Policy());
 294   return (df * scale) / (df - 2);
 295 } // mean
 296
 297 template <class RealType, class Policy>
 298 inline RealType variance(const inverse_chi_squared_distribution<RealType, Policy>& dist)
 299 { // Variance of inverse Chi-Squared distribution.
 300    RealType df = dist.degrees_of_freedom();
 301    RealType scale = dist.scale();
 302    static const char* function = "boost::math::variance(const inverse_chi_squared_distribution<%1%>&)";
 303    if(df <= 4)
 304    {
 305       return policies::raise_domain_error<RealType>(
 306          function,
 307          "inverse Chi-Squared distribution only has a variance for degrees of freedom > 4, but got degrees of freedom = %1%.",
 308          df, Policy());
 309    }
 310    return 2 * df * df * scale * scale / ((df - 2)*(df - 2) * (df - 4));
 311 } // variance
 312
 313 template <class RealType, class Policy>
 314 inline RealType mode(const inverse_chi_squared_distribution<RealType, Policy>& dist)
 315 { // mode is not defined in Mathematica.
 316   // See Discussion section http://en.wikipedia.org/wiki/Talk:Scaled-inverse-chi-square_distribution
 317   // for origin of the formula used below.
 318
 319    RealType df = dist.degrees_of_freedom();
 320    RealType scale = dist.scale();
 321    static const char* function = "boost::math::mode(const inverse_chi_squared_distribution<%1%>&)";
 322    if(df < 0)
 323       return policies::raise_domain_error<RealType>(
 324          function,
 325          "inverse Chi-Squared distribution only has a mode for degrees of freedom >= 0, but got degrees of freedom = %1%.",
 326          df, Policy());
 327    return (df * scale) / (df + 2);
 328 }
 329
 330 //template <class RealType, class Policy>
 331 //inline RealType median(const inverse_chi_squared_distribution<RealType, Policy>& dist)
 332 //{ // Median is given by Quantile[dist, 1/2]
 333 //   RealType df = dist.degrees_of_freedom();
 334 //   if(df <= 1)
 335 //      return tools::domain_error<RealType>(
 336 //         BOOST_CURRENT_FUNCTION,
 337 //         "The inverse_Chi-Squared distribution only has a median for degrees of freedom >= 0, but got degrees of freedom = %1%.",
 338 //         df);
 339 //   return df;
 340 //}
 341 // Now implemented via quantile(half) in derived accessors.
 342
 343 template <class RealType, class Policy>
 344 inline RealType skewness(const inverse_chi_squared_distribution<RealType, Policy>& dist)
 345 {
 346    BOOST_MATH_STD_USING // For ADL
 347    RealType df = dist.degrees_of_freedom();
 348    static const char* function = "boost::math::skewness(const inverse_chi_squared_distribution<%1%>&)";
 349    if(df <= 6)
 350       return policies::raise_domain_error<RealType>(
 351          function,
 352          "inverse Chi-Squared distribution only has a skewness for degrees of freedom > 6, but got degrees of freedom = %1%.",
 353          df, Policy());
 354
 355    return 4 * sqrt (2 * (df - 4)) / (df - 6);  // Not a function of scale.
 356 }
 357
 358 template <class RealType, class Policy>
 359 inline RealType kurtosis(const inverse_chi_squared_distribution<RealType, Policy>& dist)
 360 {
 361    RealType df = dist.degrees_of_freedom();
 362    static const char* function = "boost::math::kurtosis(const inverse_chi_squared_distribution<%1%>&)";
 363    if(df <= 8)
 364       return policies::raise_domain_error<RealType>(
 365          function,
 366          "inverse Chi-Squared distribution only has a kurtosis for degrees of freedom > 8, but got degrees of freedom = %1%.",
 367          df, Policy());
 368
 369    return kurtosis_excess(dist) + 3;
 370 }
 371
 372 template <class RealType, class Policy>
 373 inline RealType kurtosis_excess(const inverse_chi_squared_distribution<RealType, Policy>& dist)
 374 {
 375    RealType df = dist.degrees_of_freedom();
 376    static const char* function = "boost::math::kurtosis(const inverse_chi_squared_distribution<%1%>&)";
 377    if(df <= 8)
 378       return policies::raise_domain_error<RealType>(
 379          function,
 380          "inverse Chi-Squared distribution only has a kurtosis excess for degrees of freedom > 8, but got degrees of freedom = %1%.",
 381          df, Policy());
 382
 383    return 12 * (5 * df - 22) / ((df - 6 )*(df - 8));  // Not a function of scale.
 384 }
 385
 386 //
 387 // Parameter estimation comes last:
 388 //
 389
 390 } // namespace math
 391 } // namespace boost
 392
 393 // This include must be at the end, *after* the accessors
 394 // for this distribution have been defined, in order to
 395 // keep compilers that support two-phase lookup happy.
 396 #include <boost/math/distributions/detail/derived_accessors.hpp>
 397
 398 #endif // BOOST_MATH_DISTRIBUTIONS_INVERSE_CHI_SQUARED_HPP