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 template <class RealType, class Policy>
 101 inline const std::pair<RealType, RealType> range(const inverse_chi_squared_distribution<RealType, Policy>& /*dist*/)
 102 {  // Range of permissible values for random variable x.
 103    using boost::math::tools::max_value;
 104    return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // 0 to + infinity.
 105 }
 106
 107 template <class RealType, class Policy>
 108 inline const std::pair<RealType, RealType> support(const inverse_chi_squared_distribution<RealType, Policy>& /*dist*/)
 109 {  // Range of supported values for random variable x.
 110    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
 111    return std::pair<RealType, RealType>(static_cast<RealType>(0), tools::max_value<RealType>()); // 0 to + infinity.
 112 }
 113
 114 template <class RealType, class Policy>
 115 RealType pdf(const inverse_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
 116 {
 117    BOOST_MATH_STD_USING  // for ADL of std functions.
 118    RealType df = dist.degrees_of_freedom();
 119    RealType scale = dist.scale();
 120    RealType error_result;
 121
 122    static const char* function = "boost::math::pdf(const inverse_chi_squared_distribution<%1%>&, %1%)";
 123
 124    if(false == detail::check_inverse_chi_squared
 125      (function, df, scale, &error_result, Policy())
 126      )
 127    { // Bad distribution.
 128       return error_result;
 129    }
 130    if((x < 0) || !(boost::math::isfinite)(x))
 131    { // Bad x.
 132       return policies::raise_domain_error<RealType>(
 133          function, "inverse Chi Square parameter was %1%, but must be >= 0 !", x, Policy());
 134    }
 135
 136    if(x == 0)
 137    { // Treat as special case.
 138      return 0;
 139    }
 140    // Wikipedia scaled inverse chi sq (df, scale) related to inv gamma (df/2, df * scale /2)
 141    // so use inverse gamma pdf with shape = df/2, scale df * scale /2
 142    // RealType shape = df /2; // inv_gamma shape
 143    // RealType scale = df * scale/2; // inv_gamma scale
 144    // RealType result = gamma_p_derivative(shape, scale / x, Policy()) * scale / (x * x);
 145    RealType result = df * scale/2 / x;
 146    if(result < tools::min_value<RealType>())
 147       return 0; // Random variable is near enough infinite.
 148    result = gamma_p_derivative(df/2, result, Policy()) * df * scale/2;
 149    if(result != 0) // prevent 0 / 0,  gamma_p_derivative -> 0 faster than x^2
 150       result /= (x * x);
 151    return result;
 152 } // pdf
 153
 154 template <class RealType, class Policy>
 155 inline RealType cdf(const inverse_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
 156 {
 157    static const char* function = "boost::math::cdf(const inverse_chi_squared_distribution<%1%>&, %1%)";
 158    RealType df = dist.degrees_of_freedom();
 159    RealType scale = dist.scale();
 160    RealType error_result;
 161
 162    if(false ==
 163        detail::check_inverse_chi_squared(function, df, scale, &error_result, Policy())
 164      )
 165    { // Bad distribution.
 166       return error_result;
 167    }
 168    if((x < 0) || !(boost::math::isfinite)(x))
 169    { // Bad x.
 170       return policies::raise_domain_error<RealType>(
 171          function, "inverse Chi Square parameter was %1%, but must be >= 0 !", x, Policy());
 172    }
 173    if (x == 0)
 174    { // Treat zero as a special case.
 175      return 0;
 176    }
 177    // RealType shape = df /2; // inv_gamma shape,
 178    // RealType scale = df * scale/2; // inv_gamma scale,
 179    // result = boost::math::gamma_q(shape, scale / x, Policy()); // inverse_gamma code.
 180    return boost::math::gamma_q(df / 2, (df * (scale / 2)) / x, Policy());
 181 } // cdf
 182
 183 template <class RealType, class Policy>
 184 inline RealType quantile(const inverse_chi_squared_distribution<RealType, Policy>& dist, const RealType& p)
 185 {
 186    using boost::math::gamma_q_inv;
 187    RealType df = dist.degrees_of_freedom();
 188    RealType scale = dist.scale();
 189
 190    static const char* function = "boost::math::quantile(const inverse_chi_squared_distribution<%1%>&, %1%)";
 191    // Error check:
 192    RealType error_result;
 193    if(false == detail::check_df(
 194          function, df, &error_result, Policy())
 195          && detail::check_probability(
 196             function, p, &error_result, Policy()))
 197    {
 198       return error_result;
 199    }
 200    if(false == detail::check_probability(
 201             function, p, &error_result, Policy()))
 202    {
 203       return error_result;
 204    }
 205    // RealType shape = df /2; // inv_gamma shape,
 206    // RealType scale = df * scale/2; // inv_gamma scale,
 207    // result = scale / gamma_q_inv(shape, p, Policy());
 208       RealType result = gamma_q_inv(df /2, p, Policy());
 209       if(result == 0)
 210          return policies::raise_overflow_error<RealType, Policy>(function, "Random variable is infinite.", Policy());
 211       result = df * (scale / 2) / result;
 212       return result;
 213 } // quantile
 214
 215 template <class RealType, class Policy>
 216 inline RealType cdf(const complemented2_type<inverse_chi_squared_distribution<RealType, Policy>, RealType>& c)
 217 {
 218    using boost::math::gamma_q_inv;
 219    RealType const& df = c.dist.degrees_of_freedom();
 220    RealType const& scale = c.dist.scale();
 221    RealType const& x = c.param;
 222    static const char* function = "boost::math::cdf(const inverse_chi_squared_distribution<%1%>&, %1%)";
 223    // Error check:
 224    RealType error_result;
 225    if(false == detail::check_df(
 226          function, df, &error_result, Policy()))
 227    {
 228       return error_result;
 229    }
 230    if (x == 0)
 231    { // Treat zero as a special case.
 232      return 1;
 233    }
 234    if((x < 0) || !(boost::math::isfinite)(x))
 235    {
 236       return policies::raise_domain_error<RealType>(
 237          function, "inverse Chi Square parameter was %1%, but must be > 0 !", x, Policy());
 238    }
 239    // RealType shape = df /2; // inv_gamma shape,
 240    // RealType scale = df * scale/2; // inv_gamma scale,
 241    // result = gamma_p(shape, scale/c.param, Policy()); use inv_gamma.
 242
 243    return gamma_p(df / 2, (df * scale/2) / x, Policy()); // OK
 244 } // cdf(complemented
 245
 246 template <class RealType, class Policy>
 247 inline RealType quantile(const complemented2_type<inverse_chi_squared_distribution<RealType, Policy>, RealType>& c)
 248 {
 249    using boost::math::gamma_q_inv;
 250
 251    RealType const& df = c.dist.degrees_of_freedom();
 252    RealType const& scale = c.dist.scale();
 253    RealType const& q = c.param;
 254    static const char* function = "boost::math::quantile(const inverse_chi_squared_distribution<%1%>&, %1%)";
 255    // Error check:
 256    RealType error_result;
 257    if(false == detail::check_df(function, df, &error_result, Policy()))
 258    {
 259       return error_result;
 260    }
 261    if(false == detail::check_probability(function, q, &error_result, Policy()))
 262    {
 263       return error_result;
 264    }
 265    // RealType shape = df /2; // inv_gamma shape,
 266    // RealType scale = df * scale/2; // inv_gamma scale,
 267    // result = scale / gamma_p_inv(shape, q, Policy());  // using inv_gamma.
 268    RealType result = gamma_p_inv(df/2, q, Policy());
 269    if(result == 0)
 270       return policies::raise_overflow_error<RealType, Policy>(function, "Random variable is infinite.", Policy());
 271    result = (df * scale / 2) / result;
 272    return result;
 273 } // quantile(const complement
 274
 275 template <class RealType, class Policy>
 276 inline RealType mean(const inverse_chi_squared_distribution<RealType, Policy>& dist)
 277 { // Mean of inverse Chi-Squared distribution.
 278    RealType df = dist.degrees_of_freedom();
 279    RealType scale = dist.scale();
 280
 281    static const char* function = "boost::math::mean(const inverse_chi_squared_distribution<%1%>&)";
 282    if(df <= 2)
 283       return policies::raise_domain_error<RealType>(
 284          function,
 285          "inverse Chi-Squared distribution only has a mode for degrees of freedom > 2, but got degrees of freedom = %1%.",
 286          df, Policy());
 287   return (df * scale) / (df - 2);
 288 } // mean
 289
 290 template <class RealType, class Policy>
 291 inline RealType variance(const inverse_chi_squared_distribution<RealType, Policy>& dist)
 292 { // Variance of inverse Chi-Squared distribution.
 293    RealType df = dist.degrees_of_freedom();
 294    RealType scale = dist.scale();
 295    static const char* function = "boost::math::variance(const inverse_chi_squared_distribution<%1%>&)";
 296    if(df <= 4)
 297    {
 298       return policies::raise_domain_error<RealType>(
 299          function,
 300          "inverse Chi-Squared distribution only has a variance for degrees of freedom > 4, but got degrees of freedom = %1%.",
 301          df, Policy());
 302    }
 303    return 2 * df * df * scale * scale / ((df - 2)*(df - 2) * (df - 4));
 304 } // variance
 305
 306 template <class RealType, class Policy>
 307 inline RealType mode(const inverse_chi_squared_distribution<RealType, Policy>& dist)
 308 { // mode is not defined in Mathematica.
 309   // See Discussion section http://en.wikipedia.org/wiki/Talk:Scaled-inverse-chi-square_distribution
 310   // for origin of the formula used below.
 311
 312    RealType df = dist.degrees_of_freedom();
 313    RealType scale = dist.scale();
 314    static const char* function = "boost::math::mode(const inverse_chi_squared_distribution<%1%>&)";
 315    if(df < 0)
 316       return policies::raise_domain_error<RealType>(
 317          function,
 318          "inverse Chi-Squared distribution only has a mode for degrees of freedom >= 0, but got degrees of freedom = %1%.",
 319          df, Policy());
 320    return (df * scale) / (df + 2);
 321 }
 322
 323 //template <class RealType, class Policy>
 324 //inline RealType median(const inverse_chi_squared_distribution<RealType, Policy>& dist)
 325 //{ // Median is given by Quantile[dist, 1/2]
 326 //   RealType df = dist.degrees_of_freedom();
 327 //   if(df <= 1)
 328 //      return tools::domain_error<RealType>(
 329 //         BOOST_CURRENT_FUNCTION,
 330 //         "The inverse_Chi-Squared distribution only has a median for degrees of freedom >= 0, but got degrees of freedom = %1%.",
 331 //         df);
 332 //   return df;
 333 //}
 334 // Now implemented via quantile(half) in derived accessors.
 335
 336 template <class RealType, class Policy>
 337 inline RealType skewness(const inverse_chi_squared_distribution<RealType, Policy>& dist)
 338 {
 339    BOOST_MATH_STD_USING // For ADL
 340    RealType df = dist.degrees_of_freedom();
 341    static const char* function = "boost::math::skewness(const inverse_chi_squared_distribution<%1%>&)";
 342    if(df <= 6)
 343       return policies::raise_domain_error<RealType>(
 344          function,
 345          "inverse Chi-Squared distribution only has a skewness for degrees of freedom > 6, but got degrees of freedom = %1%.",
 346          df, Policy());
 347
 348    return 4 * sqrt (2 * (df - 4)) / (df - 6);  // Not a function of scale.
 349 }
 350
 351 template <class RealType, class Policy>
 352 inline RealType kurtosis(const inverse_chi_squared_distribution<RealType, Policy>& dist)
 353 {
 354    RealType df = dist.degrees_of_freedom();
 355    static const char* function = "boost::math::kurtosis(const inverse_chi_squared_distribution<%1%>&)";
 356    if(df <= 8)
 357       return policies::raise_domain_error<RealType>(
 358          function,
 359          "inverse Chi-Squared distribution only has a kurtosis for degrees of freedom > 8, but got degrees of freedom = %1%.",
 360          df, Policy());
 361
 362    return kurtosis_excess(dist) + 3;
 363 }
 364
 365 template <class RealType, class Policy>
 366 inline RealType kurtosis_excess(const inverse_chi_squared_distribution<RealType, Policy>& dist)
 367 {
 368    RealType df = dist.degrees_of_freedom();
 369    static const char* function = "boost::math::kurtosis(const inverse_chi_squared_distribution<%1%>&)";
 370    if(df <= 8)
 371       return policies::raise_domain_error<RealType>(
 372          function,
 373          "inverse Chi-Squared distribution only has a kurtosis excess for degrees of freedom > 8, but got degrees of freedom = %1%.",
 374          df, Policy());
 375
 376    return 12 * (5 * df - 22) / ((df - 6 )*(df - 8));  // Not a function of scale.
 377 }
 378
 379 //
 380 // Parameter estimation comes last:
 381 //
 382
 383 } // namespace math
 384 } // namespace boost
 385
 386 // This include must be at the end, *after* the accessors
 387 // for this distribution have been defined, in order to
 388 // keep compilers that support two-phase lookup happy.
 389 #include <boost/math/distributions/detail/derived_accessors.hpp>
 390
 391 #endif // BOOST_MATH_DISTRIBUTIONS_INVERSE_CHI_SQUARED_HPP