math/distributions/inverse_gamma.hpp

0001 // inverse_gamma.hpp
0002
0003 //  Copyright Paul A. Bristow 2010.
0004 //  Copyright John Maddock 2010.
0005 //  Use, modification and distribution are subject to the
0006 //  Boost Software License, Version 1.0. (See accompanying file
0007 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0008
0009 #ifndef BOOST_STATS_INVERSE_GAMMA_HPP
0010 #define BOOST_STATS_INVERSE_GAMMA_HPP
0011
0012 // Inverse Gamma Distribution is a two-parameter family
0013 // of continuous probability distributions
0014 // on the positive real line, which is the distribution of
0015 // the reciprocal of a variable distributed according to the gamma distribution.
0016
0017 // http://en.wikipedia.org/wiki/Inverse-gamma_distribution
0018 // http://rss.acs.unt.edu/Rdoc/library/pscl/html/igamma.html
0019
0020 // See also gamma distribution at gamma.hpp:
0021 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
0022 // http://mathworld.wolfram.com/GammaDistribution.html
0023 // http://en.wikipedia.org/wiki/Gamma_distribution
0024
0025 #include <boost/math/distributions/fwd.hpp>
0026 #include <boost/math/special_functions/gamma.hpp>
0027 #include <boost/math/distributions/detail/common_error_handling.hpp>
0028 #include <boost/math/distributions/complement.hpp>
0029
0030 #include <utility>
0031 #include <cfenv>
0032
0033 namespace boost{ namespace math
0034 {
0035 namespace detail
0036 {
0037
0038 template <class RealType, class Policy>
0039 inline bool check_inverse_gamma_shape(
0040       const char* function, // inverse_gamma
0041       RealType shape, // shape aka alpha
0042       RealType* result, // to update, perhaps with NaN
0043       const Policy& pol)
0044 {  // Sources say shape argument must be > 0
0045    // but seems logical to allow shape zero as special case,
0046    // returning pdf and cdf zero (but not < 0).
0047    // (Functions like mean, variance with other limits on shape are checked
0048    // in version including an operator & limit below).
0049    if((shape < 0) || !(boost::math::isfinite)(shape))
0050    {
0051       *result = policies::raise_domain_error<RealType>(
0052          function,
0053          "Shape parameter is %1%, but must be >= 0 !", shape, pol);
0054       return false;
0055    }
0056    return true;
0057 } //bool check_inverse_gamma_shape
0058
0059 template <class RealType, class Policy>
0060 inline bool check_inverse_gamma_x(
0061       const char* function,
0062       RealType const& x,
0063       RealType* result, const Policy& pol)
0064 {
0065    if((x < 0) || !(boost::math::isfinite)(x))
0066    {
0067       *result = policies::raise_domain_error<RealType>(
0068          function,
0069          "Random variate is %1% but must be >= 0 !", x, pol);
0070       return false;
0071    }
0072    return true;
0073 }
0074
0075 template <class RealType, class Policy>
0076 inline bool check_inverse_gamma(
0077       const char* function, // TODO swap these over, so shape is first.
0078       RealType scale,  // scale aka beta
0079       RealType shape, // shape aka alpha
0080       RealType* result, const Policy& pol)
0081 {
0082    return check_scale(function, scale, result, pol)
0083      && check_inverse_gamma_shape(function, shape, result, pol);
0084 } // bool check_inverse_gamma
0085
0086 } // namespace detail
0087
0088 template <class RealType = double, class Policy = policies::policy<> >
0089 class inverse_gamma_distribution
0090 {
0091 public:
0092    using value_type = RealType;
0093    using policy_type = Policy;
0094
0095    explicit inverse_gamma_distribution(RealType l_shape = 1, RealType l_scale = 1)
0096       : m_shape(l_shape), m_scale(l_scale)
0097    {
0098       RealType result;
0099       detail::check_inverse_gamma(
0100         "boost::math::inverse_gamma_distribution<%1%>::inverse_gamma_distribution",
0101         l_scale, l_shape, &result, Policy());
0102    }
0103
0104    RealType shape()const
0105    {
0106       return m_shape;
0107    }
0108
0109    RealType scale()const
0110    {
0111       return m_scale;
0112    }
0113 private:
0114    //
0115    // Data members:
0116    //
0117    RealType m_shape;     // distribution shape
0118    RealType m_scale;     // distribution scale
0119 };
0120
0121 using inverse_gamma = inverse_gamma_distribution<double>;
0122 // typedef - but potential clash with name of inverse gamma *function*.
0123 // but there is a typedef for the gamma distribution (gamma)
0124
0125 #ifdef __cpp_deduction_guides
0126 template <class RealType>
0127 inverse_gamma_distribution(RealType)->inverse_gamma_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0128 template <class RealType>
0129 inverse_gamma_distribution(RealType,RealType)->inverse_gamma_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0130 #endif
0131
0132 // Allow random variable x to be zero, treated as a special case (unlike some definitions).
0133
0134 template <class RealType, class Policy>
0135 inline std::pair<RealType, RealType> range(const inverse_gamma_distribution<RealType, Policy>& /* dist */)
0136 {  // Range of permissible values for random variable x.
0137    using boost::math::tools::max_value;
0138    return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
0139 }
0140
0141 template <class RealType, class Policy>
0142 inline std::pair<RealType, RealType> support(const inverse_gamma_distribution<RealType, Policy>& /* dist */)
0143 {  // Range of supported values for random variable x.
0144    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0145    using boost::math::tools::max_value;
0146    using boost::math::tools::min_value;
0147    return std::pair<RealType, RealType>(static_cast<RealType>(0),  max_value<RealType>());
0148 }
0149
0150 template <class RealType, class Policy>
0151 inline RealType pdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)
0152 {
0153    BOOST_MATH_STD_USING  // for ADL of std functions
0154
0155    static const char* function = "boost::math::pdf(const inverse_gamma_distribution<%1%>&, %1%)";
0156
0157    RealType shape = dist.shape();
0158    RealType scale = dist.scale();
0159
0160    RealType result = 0;
0161    if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
0162    { // distribution parameters bad.
0163       return result;
0164    }
0165    if(x == 0)
0166    { // Treat random variate zero as a special case.
0167       return 0;
0168    }
0169    else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy()))
0170    { // x bad.
0171       return result;
0172    }
0173    result = scale / x;
0174    if(result < tools::min_value<RealType>())
0175       return 0;  // random variable is infinite or so close as to make no difference.
0176    result = gamma_p_derivative(shape, result, Policy()) * scale;
0177    if(0 != result)
0178    {
0179       if(x < 0)
0180       {
0181          // x * x may under or overflow, likewise our result,
0182          // so be extra careful about the arithmetic:
0183          RealType lim = tools::max_value<RealType>() * x;
0184          if(lim < result)
0185             return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy());
0186          result /= x;
0187          if(lim < result)
0188             return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy());
0189          result /= x;
0190       }
0191       result /= (x * x);
0192    }
0193
0194    return result;
0195 } // pdf
0196
0197 template <class RealType, class Policy>
0198 inline RealType logpdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)
0199 {
0200    BOOST_MATH_STD_USING  // for ADL of std functions
0201    using boost::math::lgamma;
0202
0203    static const char* function = "boost::math::logpdf(const inverse_gamma_distribution<%1%>&, %1%)";
0204
0205    RealType shape = dist.shape();
0206    RealType scale = dist.scale();
0207
0208    RealType result = -std::numeric_limits<RealType>::infinity();
0209    if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
0210    { // distribution parameters bad.
0211       return result;
0212    }
0213    if(x == 0)
0214    { // Treat random variate zero as a special case.
0215       return result;
0216    }
0217    else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy()))
0218    { // x bad.
0219       return result;
0220    }
0221    result = scale / x;
0222    if(result < tools::min_value<RealType>())
0223       return result;  // random variable is infinite or so close as to make no difference.
0224
0225    // x * x may under or overflow, likewise our result
0226    if (!(boost::math::isfinite)(x*x))
0227    {
0228       return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy());
0229    }
0230
0231    return shape * log(scale) + (-shape-1)*log(x) - lgamma(shape) - (scale/x);
0232 } // pdf
0233
0234 template <class RealType, class Policy>
0235 inline RealType cdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)
0236 {
0237    BOOST_MATH_STD_USING  // for ADL of std functions
0238
0239    static const char* function = "boost::math::cdf(const inverse_gamma_distribution<%1%>&, %1%)";
0240
0241    RealType shape = dist.shape();
0242    RealType scale = dist.scale();
0243
0244    RealType result = 0;
0245    if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
0246    { // distribution parameters bad.
0247       return result;
0248    }
0249    if (x == 0)
0250    { // Treat zero as a special case.
0251      return 0;
0252    }
0253    else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy()))
0254    { // x bad
0255       return result;
0256    }
0257    result = boost::math::gamma_q(shape, scale / x, Policy());
0258    // result = tgamma(shape, scale / x) / tgamma(shape); // naive using tgamma
0259    return result;
0260 } // cdf
0261
0262 template <class RealType, class Policy>
0263 inline RealType quantile(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& p)
0264 {
0265    BOOST_MATH_STD_USING  // for ADL of std functions
0266    using boost::math::gamma_q_inv;
0267
0268    static const char* function = "boost::math::quantile(const inverse_gamma_distribution<%1%>&, %1%)";
0269
0270    RealType shape = dist.shape();
0271    RealType scale = dist.scale();
0272
0273    RealType result = 0;
0274    if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
0275       return result;
0276    if(false == detail::check_probability(function, p, &result, Policy()))
0277       return result;
0278    if(p == 1)
0279    {
0280       return policies::raise_overflow_error<RealType>(function, 0, Policy());
0281    }
0282    result = gamma_q_inv(shape, p, Policy());
0283    if((result < 1) && (result * tools::max_value<RealType>() < scale))
0284       return policies::raise_overflow_error<RealType, Policy>(function, "Value of random variable in inverse gamma distribution quantile is infinite.", Policy());
0285    result = scale / result;
0286    return result;
0287 }
0288
0289 template <class RealType, class Policy>
0290 inline RealType cdf(const complemented2_type<inverse_gamma_distribution<RealType, Policy>, RealType>& c)
0291 {
0292    BOOST_MATH_STD_USING  // for ADL of std functions
0293
0294    static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";
0295
0296    RealType shape = c.dist.shape();
0297    RealType scale = c.dist.scale();
0298
0299    RealType result = 0;
0300    if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
0301       return result;
0302    if(false == detail::check_inverse_gamma_x(function, c.param, &result, Policy()))
0303       return result;
0304
0305    if(c.param == 0)
0306       return 1; // Avoid division by zero
0307
0308    result = gamma_p(shape, scale/c.param, Policy());
0309    return result;
0310 }
0311
0312 template <class RealType, class Policy>
0313 inline RealType quantile(const complemented2_type<inverse_gamma_distribution<RealType, Policy>, RealType>& c)
0314 {
0315    BOOST_MATH_STD_USING  // for ADL of std functions
0316
0317    static const char* function = "boost::math::quantile(const inverse_gamma_distribution<%1%>&, %1%)";
0318
0319    RealType shape = c.dist.shape();
0320    RealType scale = c.dist.scale();
0321    RealType q = c.param;
0322
0323    RealType result = 0;
0324    if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
0325       return result;
0326    if(false == detail::check_probability(function, q, &result, Policy()))
0327       return result;
0328
0329    if(q == 0)
0330    {
0331       return policies::raise_overflow_error<RealType>(function, 0, Policy());
0332    }
0333    result = gamma_p_inv(shape, q, Policy());
0334    if((result < 1) && (result * tools::max_value<RealType>() < scale))
0335       return policies::raise_overflow_error<RealType, Policy>(function, "Value of random variable in inverse gamma distribution quantile is infinite.", Policy());
0336    result = scale / result;
0337    return result;
0338 }
0339
0340 template <class RealType, class Policy>
0341 inline RealType mean(const inverse_gamma_distribution<RealType, Policy>& dist)
0342 {
0343    BOOST_MATH_STD_USING  // for ADL of std functions
0344
0345    static const char* function = "boost::math::mean(const inverse_gamma_distribution<%1%>&)";
0346
0347    RealType shape = dist.shape();
0348    RealType scale = dist.scale();
0349
0350    RealType result = 0;
0351
0352    if(false == detail::check_scale(function, scale, &result, Policy()))
0353    {
0354      return result;
0355    }
0356    if((shape <= 1) || !(boost::math::isfinite)(shape))
0357    {
0358      result = policies::raise_domain_error<RealType>(
0359        function,
0360        "Shape parameter is %1%, but for a defined mean it must be > 1", shape, Policy());
0361      return result;
0362    }
0363   result = scale / (shape - 1);
0364   return result;
0365 } // mean
0366
0367 template <class RealType, class Policy>
0368 inline RealType variance(const inverse_gamma_distribution<RealType, Policy>& dist)
0369 {
0370    BOOST_MATH_STD_USING  // for ADL of std functions
0371
0372    static const char* function = "boost::math::variance(const inverse_gamma_distribution<%1%>&)";
0373
0374    RealType shape = dist.shape();
0375    RealType scale = dist.scale();
0376
0377    RealType result = 0;
0378       if(false == detail::check_scale(function, scale, &result, Policy()))
0379    {
0380      return result;
0381    }
0382    if((shape <= 2) || !(boost::math::isfinite)(shape))
0383    {
0384      result = policies::raise_domain_error<RealType>(
0385        function,
0386        "Shape parameter is %1%, but for a defined variance it must be > 2", shape, Policy());
0387      return result;
0388    }
0389    result = (scale * scale) / ((shape - 1) * (shape -1) * (shape -2));
0390    return result;
0391 }
0392
0393 template <class RealType, class Policy>
0394 inline RealType mode(const inverse_gamma_distribution<RealType, Policy>& dist)
0395 {
0396    BOOST_MATH_STD_USING  // for ADL of std functions
0397
0398    static const char* function = "boost::math::mode(const inverse_gamma_distribution<%1%>&)";
0399
0400    RealType shape = dist.shape();
0401    RealType scale = dist.scale();
0402
0403    RealType result = 0;
0404    if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))
0405    {
0406       return result;
0407    }
0408    // Only defined for shape >= 0, but is checked by check_inverse_gamma.
0409    result = scale / (shape + 1);
0410    return result;
0411 }
0412
0413 //template <class RealType, class Policy>
0414 //inline RealType median(const gamma_distribution<RealType, Policy>& dist)
0415 //{  // Wikipedia does not define median,
0416      // so rely on default definition quantile(0.5) in derived accessors.
0417 //  return result.
0418 //}
0419
0420 template <class RealType, class Policy>
0421 inline RealType skewness(const inverse_gamma_distribution<RealType, Policy>& dist)
0422 {
0423    BOOST_MATH_STD_USING  // for ADL of std functions
0424
0425    static const char* function = "boost::math::skewness(const inverse_gamma_distribution<%1%>&)";
0426
0427    RealType shape = dist.shape();
0428    RealType scale = dist.scale();
0429    RealType result = 0;
0430
0431    if(false == detail::check_scale(function, scale, &result, Policy()))
0432    {
0433      return result;
0434    }
0435    if((shape <= 3) || !(boost::math::isfinite)(shape))
0436    {
0437      result = policies::raise_domain_error<RealType>(
0438        function,
0439        "Shape parameter is %1%, but for a defined skewness it must be > 3", shape, Policy());
0440      return result;
0441    }
0442    result = (4 * sqrt(shape - 2) ) / (shape - 3);
0443    return result;
0444 }
0445
0446 template <class RealType, class Policy>
0447 inline RealType kurtosis_excess(const inverse_gamma_distribution<RealType, Policy>& dist)
0448 {
0449    BOOST_MATH_STD_USING  // for ADL of std functions
0450
0451    static const char* function = "boost::math::kurtosis_excess(const inverse_gamma_distribution<%1%>&)";
0452
0453    RealType shape = dist.shape();
0454    RealType scale = dist.scale();
0455
0456    RealType result = 0;
0457    if(false == detail::check_scale(function, scale, &result, Policy()))
0458    {
0459      return result;
0460    }
0461    if((shape <= 4) || !(boost::math::isfinite)(shape))
0462    {
0463      result = policies::raise_domain_error<RealType>(
0464        function,
0465        "Shape parameter is %1%, but for a defined kurtosis excess it must be > 4", shape, Policy());
0466      return result;
0467    }
0468    result = (30 * shape - 66) / ((shape - 3) * (shape - 4));
0469    return result;
0470 }
0471
0472 template <class RealType, class Policy>
0473 inline RealType kurtosis(const inverse_gamma_distribution<RealType, Policy>& dist)
0474 {
0475   static const char* function = "boost::math::kurtosis(const inverse_gamma_distribution<%1%>&)";
0476    RealType shape = dist.shape();
0477    RealType scale = dist.scale();
0478
0479    RealType result = 0;
0480
0481   if(false == detail::check_scale(function, scale, &result, Policy()))
0482    {
0483      return result;
0484    }
0485    if((shape <= 4) || !(boost::math::isfinite)(shape))
0486    {
0487      result = policies::raise_domain_error<RealType>(
0488        function,
0489        "Shape parameter is %1%, but for a defined kurtosis it must be > 4", shape, Policy());
0490      return result;
0491    }
0492   return kurtosis_excess(dist) + 3;
0493 }
0494
0495 } // namespace math
0496 } // namespace boost
0497
0498 // This include must be at the end, *after* the accessors
0499 // for this distribution have been defined, in order to
0500 // keep compilers that support two-phase lookup happy.
0501 #include <boost/math/distributions/detail/derived_accessors.hpp>
0502
0503 #endif // BOOST_STATS_INVERSE_GAMMA_HPP