ceph/src/boost/libs/math/include/boost/math/distributions/inverse_gamma.hpp

   1 // inverse_gamma.hpp
   2
   3 //  Copyright Paul A. Bristow 2010.
   4 //  Copyright John Maddock 2010.
   5 //  Use, modification and distribution are subject to the
   6 //  Boost Software License, Version 1.0. (See accompanying file
   7 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   8
   9 #ifndef BOOST_STATS_INVERSE_GAMMA_HPP
  10 #define BOOST_STATS_INVERSE_GAMMA_HPP
  11
  12 // Inverse Gamma Distribution is a two-parameter family
  13 // of continuous probability distributions
  14 // on the positive real line, which is the distribution of
  15 // the reciprocal of a variable distributed according to the gamma distribution.
  16
  17 // http://en.wikipedia.org/wiki/Inverse-gamma_distribution
  18 // http://rss.acs.unt.edu/Rdoc/library/pscl/html/igamma.html
  19
  20 // See also gamma distribution at gamma.hpp:
  21 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
  22 // http://mathworld.wolfram.com/GammaDistribution.html
  23 // http://en.wikipedia.org/wiki/Gamma_distribution
  24
  25 #include <boost/math/distributions/fwd.hpp>
  26 #include <boost/math/special_functions/gamma.hpp>
  27 #include <boost/math/distributions/detail/common_error_handling.hpp>
  28 #include <boost/math/distributions/complement.hpp>
  29
  30 #include <utility>
  31
  32 namespace boost{ namespace math
  33 {
  34 namespace detail
  35 {
  36
  37 template <class RealType, class Policy>
  38 inline bool check_inverse_gamma_shape(
  39       const char* function, // inverse_gamma
  40       RealType shape, // shape aka alpha
  41       RealType* result, // to update, perhaps with NaN
  42       const Policy& pol)
  43 {  // Sources say shape argument must be > 0
  44    // but seems logical to allow shape zero as special case,
  45    // returning pdf and cdf zero (but not < 0).
  46    // (Functions like mean, variance with other limits on shape are checked
  47    // in version including an operator & limit below).
  48    if((shape < 0) || !(boost::math::isfinite)(shape))
  49    {
  50       *result = policies::raise_domain_error<RealType>(
  51          function,
  52          "Shape parameter is %1%, but must be >= 0 !", shape, pol);
  53       return false;
  54    }
  55    return true;
  56 } //bool check_inverse_gamma_shape
  57
  58 template <class RealType, class Policy>
  59 inline bool check_inverse_gamma_x(
  60       const char* function,
  61       RealType const& x,
  62       RealType* result, const Policy& pol)
  63 {
  64    if((x < 0) || !(boost::math::isfinite)(x))
  65    {
  66       *result = policies::raise_domain_error<RealType>(
  67          function,
  68          "Random variate is %1% but must be >= 0 !", x, pol);
  69       return false;
  70    }
  71    return true;
  72 }
  73
  74 template <class RealType, class Policy>
  75 inline bool check_inverse_gamma(
  76       const char* function, // TODO swap these over, so shape is first.
  77       RealType scale,  // scale aka beta
  78       RealType shape, // shape aka alpha
  79       RealType* result, const Policy& pol)
  80 {
  81    return check_scale(function, scale, result, pol)
  82      && check_inverse_gamma_shape(function, shape, result, pol);
  83 } // bool check_inverse_gamma
  84
  85 } // namespace detail
  86
  87 template <class RealType = double, class Policy = policies::policy<> >
  88 class inverse_gamma_distribution
  89 {
  90 public:
  91    typedef RealType value_type;
  92    typedef Policy policy_type;
  93
  94    inverse_gamma_distribution(RealType l_shape = 1, RealType l_scale = 1)
  95       : m_shape(l_shape), m_scale(l_scale)
  96    {
  97       RealType result;
  98       detail::check_inverse_gamma(
  99         "boost::math::inverse_gamma_distribution<%1%>::inverse_gamma_distribution",
 100         l_scale, l_shape, &result, Policy());
 101    }
 102
 103    RealType shape()const
 104    {
 105       return m_shape;
 106    }
 107
 108    RealType scale()const
 109    {
 110       return m_scale;
 111    }
 112 private:
 113    //
 114    // Data members:
 115    //
 116    RealType m_shape;     // distribution shape
 117    RealType m_scale;     // distribution scale
 118 };
 119
 120 typedef inverse_gamma_distribution<double> inverse_gamma;
 121 // typedef - but potential clash with name of inverse gamma *function*.
 122 // but there is a typedef for gamma
 123 //   typedef boost::math::gamma_distribution<Type, Policy> gamma;
 124
 125 // Allow random variable x to be zero, treated as a special case (unlike some definitions).
 126
 127 template <class RealType, class Policy>
 128 inline const std::pair<RealType, RealType> range(const inverse_gamma_distribution<RealType, Policy>& /* dist */)
 129 {  // Range of permissible values for random variable x.
 130    using boost::math::tools::max_value;
 131    return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
 132 }
 133
 134 template <class RealType, class Policy>
 135 inline const std::pair<RealType, RealType> support(const inverse_gamma_distribution<RealType, Policy>& /* dist */)
 136 {  // Range of supported values for random variable x.
 137    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
 138    using boost::math::tools::max_value;
 139    using boost::math::tools::min_value;
 140    return std::pair<RealType, RealType>(static_cast<RealType>(0),  max_value<RealType>());
 141 }
 142
 143 template <class RealType, class Policy>
 144 inline RealType pdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)
 145 {
 146    BOOST_MATH_STD_USING  // for ADL of std functions
 147
 148    static const char* function = "boost::math::pdf(const inverse_gamma_distribution<%1%>&, %1%)";
 149
 150    RealType shape = dist.shape();
 151    RealType scale = dist.scale();
 152
 153    RealType result = 0;
 154    if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
 155    { // distribution parameters bad.
 156       return result;
 157    }
 158    if(x == 0)
 159    { // Treat random variate zero as a special case.
 160       return 0;
 161    }
 162    else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy()))
 163    { // x bad.
 164       return result;
 165    }
 166    result = scale / x;
 167    if(result < tools::min_value<RealType>())
 168       return 0;  // random variable is infinite or so close as to make no difference.
 169    result = gamma_p_derivative(shape, result, Policy()) * scale;
 170    if(0 != result)
 171    {
 172       if(x < 0)
 173       {
 174          // x * x may under or overflow, likewise our result,
 175          // so be extra careful about the arithmetic:
 176          RealType lim = tools::max_value<RealType>() * x;
 177          if(lim < result)
 178             return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy());
 179          result /= x;
 180          if(lim < result)
 181             return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy());
 182          result /= x;
 183       }
 184       result /= (x * x);
 185    }
 186    // better than naive
 187    // result = (pow(scale, shape) * pow(x, (-shape -1)) * exp(-scale/x) ) / tgamma(shape);
 188    return result;
 189 } // pdf
 190
 191 template <class RealType, class Policy>
 192 inline RealType cdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)
 193 {
 194    BOOST_MATH_STD_USING  // for ADL of std functions
 195
 196    static const char* function = "boost::math::cdf(const inverse_gamma_distribution<%1%>&, %1%)";
 197
 198    RealType shape = dist.shape();
 199    RealType scale = dist.scale();
 200
 201    RealType result = 0;
 202    if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
 203    { // distribution parameters bad.
 204       return result;
 205    }
 206    if (x == 0)
 207    { // Treat zero as a special case.
 208      return 0;
 209    }
 210    else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy()))
 211    { // x bad
 212       return result;
 213    }
 214    result = boost::math::gamma_q(shape, scale / x, Policy());
 215    // result = tgamma(shape, scale / x) / tgamma(shape); // naive using tgamma
 216    return result;
 217 } // cdf
 218
 219 template <class RealType, class Policy>
 220 inline RealType quantile(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& p)
 221 {
 222    BOOST_MATH_STD_USING  // for ADL of std functions
 223    using boost::math::gamma_q_inv;
 224
 225    static const char* function = "boost::math::quantile(const inverse_gamma_distribution<%1%>&, %1%)";
 226
 227    RealType shape = dist.shape();
 228    RealType scale = dist.scale();
 229
 230    RealType result = 0;
 231    if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
 232       return result;
 233    if(false == detail::check_probability(function, p, &result, Policy()))
 234       return result;
 235    if(p == 1)
 236    {
 237       return policies::raise_overflow_error<RealType>(function, 0, Policy());
 238    }
 239    result = gamma_q_inv(shape, p, Policy());
 240    if((result < 1) && (result * tools::max_value<RealType>() < scale))
 241       return policies::raise_overflow_error<RealType, Policy>(function, "Value of random variable in inverse gamma distribution quantile is infinite.", Policy());
 242    result = scale / result;
 243    return result;
 244 }
 245
 246 template <class RealType, class Policy>
 247 inline RealType cdf(const complemented2_type<inverse_gamma_distribution<RealType, Policy>, RealType>& c)
 248 {
 249    BOOST_MATH_STD_USING  // for ADL of std functions
 250
 251    static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";
 252
 253    RealType shape = c.dist.shape();
 254    RealType scale = c.dist.scale();
 255
 256    RealType result = 0;
 257    if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
 258       return result;
 259    if(false == detail::check_inverse_gamma_x(function, c.param, &result, Policy()))
 260       return result;
 261
 262    if(c.param == 0)
 263       return 1; // Avoid division by zero
 264
 265    //result = 1. - gamma_q(shape, c.param / scale, Policy());
 266    result = gamma_p(shape, scale/c.param, Policy());
 267    return result;
 268 }
 269
 270 template <class RealType, class Policy>
 271 inline RealType quantile(const complemented2_type<inverse_gamma_distribution<RealType, Policy>, RealType>& c)
 272 {
 273    BOOST_MATH_STD_USING  // for ADL of std functions
 274
 275    static const char* function = "boost::math::quantile(const inverse_gamma_distribution<%1%>&, %1%)";
 276
 277    RealType shape = c.dist.shape();
 278    RealType scale = c.dist.scale();
 279    RealType q = c.param;
 280
 281    RealType result = 0;
 282    if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
 283       return result;
 284    if(false == detail::check_probability(function, q, &result, Policy()))
 285       return result;
 286
 287    if(q == 0)
 288    {
 289       return policies::raise_overflow_error<RealType>(function, 0, Policy());
 290    }
 291    result = gamma_p_inv(shape, q, Policy());
 292    if((result < 1) && (result * tools::max_value<RealType>() < scale))
 293       return policies::raise_overflow_error<RealType, Policy>(function, "Value of random variable in inverse gamma distribution quantile is infinite.", Policy());
 294    result = scale / result;
 295    return result;
 296 }
 297
 298 template <class RealType, class Policy>
 299 inline RealType mean(const inverse_gamma_distribution<RealType, Policy>& dist)
 300 {
 301    BOOST_MATH_STD_USING  // for ADL of std functions
 302
 303    static const char* function = "boost::math::mean(const inverse_gamma_distribution<%1%>&)";
 304
 305    RealType shape = dist.shape();
 306    RealType scale = dist.scale();
 307
 308    RealType result = 0;
 309
 310    if(false == detail::check_scale(function, scale, &result, Policy()))
 311    {
 312      return result;
 313    }
 314    if((shape <= 1) || !(boost::math::isfinite)(shape))
 315    {
 316      result = policies::raise_domain_error<RealType>(
 317        function,
 318        "Shape parameter is %1%, but for a defined mean it must be > 1", shape, Policy());
 319      return result;
 320    }
 321   result = scale / (shape - 1);
 322   return result;
 323 } // mean
 324
 325 template <class RealType, class Policy>
 326 inline RealType variance(const inverse_gamma_distribution<RealType, Policy>& dist)
 327 {
 328    BOOST_MATH_STD_USING  // for ADL of std functions
 329
 330    static const char* function = "boost::math::variance(const inverse_gamma_distribution<%1%>&)";
 331
 332    RealType shape = dist.shape();
 333    RealType scale = dist.scale();
 334
 335    RealType result = 0;
 336       if(false == detail::check_scale(function, scale, &result, Policy()))
 337    {
 338      return result;
 339    }
 340    if((shape <= 2) || !(boost::math::isfinite)(shape))
 341    {
 342      result = policies::raise_domain_error<RealType>(
 343        function,
 344        "Shape parameter is %1%, but for a defined variance it must be > 2", shape, Policy());
 345      return result;
 346    }
 347    result = (scale * scale) / ((shape - 1) * (shape -1) * (shape -2));
 348    return result;
 349 }
 350
 351 template <class RealType, class Policy>
 352 inline RealType mode(const inverse_gamma_distribution<RealType, Policy>& dist)
 353 {
 354    BOOST_MATH_STD_USING  // for ADL of std functions
 355
 356    static const char* function = "boost::math::mode(const inverse_gamma_distribution<%1%>&)";
 357
 358    RealType shape = dist.shape();
 359    RealType scale = dist.scale();
 360
 361    RealType result = 0;
 362    if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
 363    {
 364       return result;
 365    }
 366    // Only defined for shape >= 0, but is checked by check_inverse_gamma.
 367    result = scale / (shape + 1);
 368    return result;
 369 }
 370
 371 //template <class RealType, class Policy>
 372 //inline RealType median(const gamma_distribution<RealType, Policy>& dist)
 373 //{  // Wikipedia does not define median,
 374      // so rely on default definition quantile(0.5) in derived accessors.
 375 //  return result.
 376 //}
 377
 378 template <class RealType, class Policy>
 379 inline RealType skewness(const inverse_gamma_distribution<RealType, Policy>& dist)
 380 {
 381    BOOST_MATH_STD_USING  // for ADL of std functions
 382
 383    static const char* function = "boost::math::skewness(const inverse_gamma_distribution<%1%>&)";
 384
 385    RealType shape = dist.shape();
 386    RealType scale = dist.scale();
 387    RealType result = 0;
 388
 389    if(false == detail::check_scale(function, scale, &result, Policy()))
 390    {
 391      return result;
 392    }
 393    if((shape <= 3) || !(boost::math::isfinite)(shape))
 394    {
 395      result = policies::raise_domain_error<RealType>(
 396        function,
 397        "Shape parameter is %1%, but for a defined skewness it must be > 3", shape, Policy());
 398      return result;
 399    }
 400    result = (4 * sqrt(shape - 2) ) / (shape - 3);
 401    return result;
 402 }
 403
 404 template <class RealType, class Policy>
 405 inline RealType kurtosis_excess(const inverse_gamma_distribution<RealType, Policy>& dist)
 406 {
 407    BOOST_MATH_STD_USING  // for ADL of std functions
 408
 409    static const char* function = "boost::math::kurtosis_excess(const inverse_gamma_distribution<%1%>&)";
 410
 411    RealType shape = dist.shape();
 412    RealType scale = dist.scale();
 413
 414    RealType result = 0;
 415    if(false == detail::check_scale(function, scale, &result, Policy()))
 416    {
 417      return result;
 418    }
 419    if((shape <= 4) || !(boost::math::isfinite)(shape))
 420    {
 421      result = policies::raise_domain_error<RealType>(
 422        function,
 423        "Shape parameter is %1%, but for a defined kurtosis excess it must be > 4", shape, Policy());
 424      return result;
 425    }
 426    result = (30 * shape - 66) / ((shape - 3) * (shape - 4));
 427    return result;
 428 }
 429
 430 template <class RealType, class Policy>
 431 inline RealType kurtosis(const inverse_gamma_distribution<RealType, Policy>& dist)
 432 {
 433   static const char* function = "boost::math::kurtosis(const inverse_gamma_distribution<%1%>&)";
 434    RealType shape = dist.shape();
 435    RealType scale = dist.scale();
 436
 437    RealType result = 0;
 438
 439   if(false == detail::check_scale(function, scale, &result, Policy()))
 440    {
 441      return result;
 442    }
 443    if((shape <= 4) || !(boost::math::isfinite)(shape))
 444    {
 445      result = policies::raise_domain_error<RealType>(
 446        function,
 447        "Shape parameter is %1%, but for a defined kurtosis it must be > 4", shape, Policy());
 448      return result;
 449    }
 450   return kurtosis_excess(dist) + 3;
 451 }
 452
 453 } // namespace math
 454 } // namespace boost
 455
 456 // This include must be at the end, *after* the accessors
 457 // for this distribution have been defined, in order to
 458 // keep compilers that support two-phase lookup happy.
 459 #include <boost/math/distributions/detail/derived_accessors.hpp>
 460
 461 #endif // BOOST_STATS_INVERSE_GAMMA_HPP