math/distributions/inverse_gamma.hpp

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