505 lines · plain
1// inverse_gamma.hpp2 3// Copyright Paul A. Bristow 2010.4// Copyright John Maddock 2010.5// Copyright Matt Borland 2024.6// Use, modification and distribution are subject to the7// Boost Software License, Version 1.0. (See accompanying file8// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)9 10#ifndef BOOST_STATS_INVERSE_GAMMA_HPP11#define BOOST_STATS_INVERSE_GAMMA_HPP12 13// Inverse Gamma Distribution is a two-parameter family14// of continuous probability distributions15// on the positive real line, which is the distribution of16// the reciprocal of a variable distributed according to the gamma distribution.17 18// http://en.wikipedia.org/wiki/Inverse-gamma_distribution19// http://rss.acs.unt.edu/Rdoc/library/pscl/html/igamma.html20 21// See also gamma distribution at gamma.hpp:22// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm23// http://mathworld.wolfram.com/GammaDistribution.html24// http://en.wikipedia.org/wiki/Gamma_distribution25 26#include <boost/math/tools/config.hpp>27#include <boost/math/tools/tuple.hpp>28#include <boost/math/tools/numeric_limits.hpp>29#include <boost/math/distributions/fwd.hpp>30#include <boost/math/special_functions/gamma.hpp>31#include <boost/math/distributions/detail/common_error_handling.hpp>32#include <boost/math/distributions/complement.hpp>33 34namespace boost{ namespace math35{36namespace detail37{38 39template <class RealType, class Policy>40BOOST_MATH_GPU_ENABLED inline bool check_inverse_gamma_shape(41 const char* function, // inverse_gamma42 RealType shape, // shape aka alpha43 RealType* result, // to update, perhaps with NaN44 const Policy& pol)45{ // Sources say shape argument must be > 046 // but seems logical to allow shape zero as special case,47 // returning pdf and cdf zero (but not < 0).48 // (Functions like mean, variance with other limits on shape are checked49 // in version including an operator & limit below).50 if((shape < 0) || !(boost::math::isfinite)(shape))51 {52 *result = policies::raise_domain_error<RealType>(53 function,54 "Shape parameter is %1%, but must be >= 0 !", shape, pol);55 return false;56 }57 return true;58} //bool check_inverse_gamma_shape59 60template <class RealType, class Policy>61BOOST_MATH_GPU_ENABLED inline bool check_inverse_gamma_x(62 const char* function,63 RealType const& x,64 RealType* result, const Policy& pol)65{66 if((x < 0) || !(boost::math::isfinite)(x))67 {68 *result = policies::raise_domain_error<RealType>(69 function,70 "Random variate is %1% but must be >= 0 !", x, pol);71 return false;72 }73 return true;74}75 76template <class RealType, class Policy>77BOOST_MATH_GPU_ENABLED inline bool check_inverse_gamma(78 const char* function, // TODO swap these over, so shape is first.79 RealType scale, // scale aka beta80 RealType shape, // shape aka alpha81 RealType* result, const Policy& pol)82{83 return check_scale(function, scale, result, pol)84 && check_inverse_gamma_shape(function, shape, result, pol);85} // bool check_inverse_gamma86 87} // namespace detail88 89template <class RealType = double, class Policy = policies::policy<> >90class inverse_gamma_distribution91{92public:93 using value_type = RealType;94 using policy_type = Policy;95 96 BOOST_MATH_GPU_ENABLED explicit inverse_gamma_distribution(RealType l_shape = 1, RealType l_scale = 1)97 : m_shape(l_shape), m_scale(l_scale)98 {99 RealType result;100 detail::check_inverse_gamma(101 "boost::math::inverse_gamma_distribution<%1%>::inverse_gamma_distribution",102 l_scale, l_shape, &result, Policy());103 }104 105 BOOST_MATH_GPU_ENABLED RealType shape()const106 {107 return m_shape;108 }109 110 BOOST_MATH_GPU_ENABLED RealType scale()const111 {112 return m_scale;113 }114private:115 //116 // Data members:117 //118 RealType m_shape; // distribution shape119 RealType m_scale; // distribution scale120};121 122using inverse_gamma = inverse_gamma_distribution<double>;123// typedef - but potential clash with name of inverse gamma *function*.124// but there is a typedef for the gamma distribution (gamma)125 126#ifdef __cpp_deduction_guides127template <class RealType>128inverse_gamma_distribution(RealType)->inverse_gamma_distribution<typename boost::math::tools::promote_args<RealType>::type>;129template <class RealType>130inverse_gamma_distribution(RealType,RealType)->inverse_gamma_distribution<typename boost::math::tools::promote_args<RealType>::type>;131#endif132 133// Allow random variable x to be zero, treated as a special case (unlike some definitions).134 135template <class RealType, class Policy>136BOOST_MATH_GPU_ENABLED inline boost::math::pair<RealType, RealType> range(const inverse_gamma_distribution<RealType, Policy>& /* dist */)137{ // Range of permissible values for random variable x.138 using boost::math::tools::max_value;139 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());140}141 142template <class RealType, class Policy>143BOOST_MATH_GPU_ENABLED inline boost::math::pair<RealType, RealType> support(const inverse_gamma_distribution<RealType, Policy>& /* dist */)144{ // Range of supported values for random variable x.145 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.146 using boost::math::tools::max_value;147 using boost::math::tools::min_value;148 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());149}150 151template <class RealType, class Policy>152BOOST_MATH_GPU_ENABLED inline RealType pdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)153{154 BOOST_MATH_STD_USING // for ADL of std functions155 156 constexpr auto function = "boost::math::pdf(const inverse_gamma_distribution<%1%>&, %1%)";157 158 RealType shape = dist.shape();159 RealType scale = dist.scale();160 161 RealType result = 0;162 if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))163 { // distribution parameters bad.164 return result;165 } 166 if(x == 0)167 { // Treat random variate zero as a special case.168 return 0;169 }170 else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy()))171 { // x bad.172 return result;173 }174 result = scale / x;175 if(result < tools::min_value<RealType>())176 return 0; // random variable is infinite or so close as to make no difference.177 result = gamma_p_derivative(shape, result, Policy()) * scale;178 if(0 != result)179 {180 if(x < 0)181 {182 // x * x may under or overflow, likewise our result,183 // so be extra careful about the arithmetic:184 RealType lim = tools::max_value<RealType>() * x;185 if(lim < result)186 return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy());187 result /= x;188 if(lim < result)189 return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy());190 result /= x;191 }192 result /= (x * x);193 }194 195 return result;196} // pdf197 198template <class RealType, class Policy>199BOOST_MATH_GPU_ENABLED inline RealType logpdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)200{201 BOOST_MATH_STD_USING // for ADL of std functions202 using boost::math::lgamma;203 204 constexpr auto function = "boost::math::logpdf(const inverse_gamma_distribution<%1%>&, %1%)";205 206 RealType shape = dist.shape();207 RealType scale = dist.scale();208 209 RealType result = -boost::math::numeric_limits<RealType>::infinity();210 if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))211 { // distribution parameters bad.212 return result;213 } 214 if(x == 0)215 { // Treat random variate zero as a special case.216 return result;217 }218 else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy()))219 { // x bad.220 return result;221 }222 result = scale / x;223 if(result < tools::min_value<RealType>())224 return result; // random variable is infinite or so close as to make no difference.225 226 // x * x may under or overflow, likewise our result227 if (!(boost::math::isfinite)(x*x))228 {229 return policies::raise_overflow_error<RealType, Policy>(function, "PDF is infinite.", Policy());230 }231 232 return shape * log(scale) + (-shape-1)*log(x) - lgamma(shape) - (scale/x);233} // pdf234 235template <class RealType, class Policy>236BOOST_MATH_GPU_ENABLED inline RealType cdf(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& x)237{238 BOOST_MATH_STD_USING // for ADL of std functions239 240 constexpr auto function = "boost::math::cdf(const inverse_gamma_distribution<%1%>&, %1%)";241 242 RealType shape = dist.shape();243 RealType scale = dist.scale();244 245 RealType result = 0;246 if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))247 { // distribution parameters bad.248 return result;249 }250 if (x == 0)251 { // Treat zero as a special case.252 return 0;253 }254 else if(false == detail::check_inverse_gamma_x(function, x, &result, Policy()))255 { // x bad256 return result;257 }258 result = boost::math::gamma_q(shape, scale / x, Policy());259 // result = tgamma(shape, scale / x) / tgamma(shape); // naive using tgamma260 return result;261} // cdf262 263template <class RealType, class Policy>264BOOST_MATH_GPU_ENABLED inline RealType quantile(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& p)265{266 BOOST_MATH_STD_USING // for ADL of std functions267 using boost::math::gamma_q_inv;268 269 constexpr auto function = "boost::math::quantile(const inverse_gamma_distribution<%1%>&, %1%)";270 271 RealType shape = dist.shape();272 RealType scale = dist.scale();273 274 RealType result = 0;275 if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))276 return result;277 if(false == detail::check_probability(function, p, &result, Policy()))278 return result;279 if(p == 1)280 {281 return policies::raise_overflow_error<RealType>(function, 0, Policy());282 }283 result = gamma_q_inv(shape, p, Policy());284 if((result < 1) && (result * tools::max_value<RealType>() < scale))285 return policies::raise_overflow_error<RealType, Policy>(function, "Value of random variable in inverse gamma distribution quantile is infinite.", Policy());286 result = scale / result;287 return result;288}289 290template <class RealType, class Policy>291BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<inverse_gamma_distribution<RealType, Policy>, RealType>& c)292{293 BOOST_MATH_STD_USING // for ADL of std functions294 295 constexpr auto function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";296 297 RealType shape = c.dist.shape();298 RealType scale = c.dist.scale();299 300 RealType result = 0;301 if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))302 return result;303 if(false == detail::check_inverse_gamma_x(function, c.param, &result, Policy()))304 return result;305 306 if(c.param == 0)307 return 1; // Avoid division by zero308 309 result = gamma_p(shape, scale/c.param, Policy());310 return result;311}312 313template <class RealType, class Policy>314BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<inverse_gamma_distribution<RealType, Policy>, RealType>& c)315{316 BOOST_MATH_STD_USING // for ADL of std functions317 318 constexpr auto function = "boost::math::quantile(const inverse_gamma_distribution<%1%>&, %1%)";319 320 RealType shape = c.dist.shape();321 RealType scale = c.dist.scale();322 RealType q = c.param;323 324 RealType result = 0;325 if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))326 return result;327 if(false == detail::check_probability(function, q, &result, Policy()))328 return result;329 330 if(q == 0)331 {332 return policies::raise_overflow_error<RealType>(function, 0, Policy());333 }334 result = gamma_p_inv(shape, q, Policy());335 if((result < 1) && (result * tools::max_value<RealType>() < scale))336 return policies::raise_overflow_error<RealType, Policy>(function, "Value of random variable in inverse gamma distribution quantile is infinite.", Policy());337 result = scale / result;338 return result;339}340 341template <class RealType, class Policy>342BOOST_MATH_GPU_ENABLED inline RealType mean(const inverse_gamma_distribution<RealType, Policy>& dist)343{344 BOOST_MATH_STD_USING // for ADL of std functions345 346 constexpr auto function = "boost::math::mean(const inverse_gamma_distribution<%1%>&)";347 348 RealType shape = dist.shape();349 RealType scale = dist.scale();350 351 RealType result = 0;352 353 if(false == detail::check_scale(function, scale, &result, Policy()))354 {355 return result;356 }357 if((shape <= 1) || !(boost::math::isfinite)(shape))358 {359 result = policies::raise_domain_error<RealType>(360 function,361 "Shape parameter is %1%, but for a defined mean it must be > 1", shape, Policy());362 return result;363 }364 result = scale / (shape - 1);365 return result;366} // mean367 368template <class RealType, class Policy>369BOOST_MATH_GPU_ENABLED inline RealType variance(const inverse_gamma_distribution<RealType, Policy>& dist)370{371 BOOST_MATH_STD_USING // for ADL of std functions372 373 constexpr auto function = "boost::math::variance(const inverse_gamma_distribution<%1%>&)";374 375 RealType shape = dist.shape();376 RealType scale = dist.scale();377 378 RealType result = 0;379 if(false == detail::check_scale(function, scale, &result, Policy()))380 {381 return result;382 }383 if((shape <= 2) || !(boost::math::isfinite)(shape))384 {385 result = policies::raise_domain_error<RealType>(386 function,387 "Shape parameter is %1%, but for a defined variance it must be > 2", shape, Policy());388 return result;389 }390 result = (scale * scale) / ((shape - 1) * (shape -1) * (shape -2));391 return result;392}393 394template <class RealType, class Policy>395BOOST_MATH_GPU_ENABLED inline RealType mode(const inverse_gamma_distribution<RealType, Policy>& dist)396{397 BOOST_MATH_STD_USING // for ADL of std functions398 399 constexpr auto function = "boost::math::mode(const inverse_gamma_distribution<%1%>&)";400 401 RealType shape = dist.shape();402 RealType scale = dist.scale();403 404 RealType result = 0;405 if(false == detail::check_inverse_gamma(function, scale, shape, &result, Policy()))406 {407 return result;408 }409 // Only defined for shape >= 0, but is checked by check_inverse_gamma.410 result = scale / (shape + 1);411 return result;412}413 414//template <class RealType, class Policy>415//inline RealType median(const gamma_distribution<RealType, Policy>& dist)416//{ // Wikipedia does not define median,417 // so rely on default definition quantile(0.5) in derived accessors.418// return result.419//}420 421template <class RealType, class Policy>422BOOST_MATH_GPU_ENABLED inline RealType skewness(const inverse_gamma_distribution<RealType, Policy>& dist)423{424 BOOST_MATH_STD_USING // for ADL of std functions425 426 constexpr auto function = "boost::math::skewness(const inverse_gamma_distribution<%1%>&)";427 428 RealType shape = dist.shape();429 RealType scale = dist.scale();430 RealType result = 0;431 432 if(false == detail::check_scale(function, scale, &result, Policy()))433 {434 return result;435 }436 if((shape <= 3) || !(boost::math::isfinite)(shape))437 {438 result = policies::raise_domain_error<RealType>(439 function,440 "Shape parameter is %1%, but for a defined skewness it must be > 3", shape, Policy());441 return result;442 }443 result = (4 * sqrt(shape - 2) ) / (shape - 3);444 return result;445}446 447template <class RealType, class Policy>448BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const inverse_gamma_distribution<RealType, Policy>& dist)449{450 BOOST_MATH_STD_USING // for ADL of std functions451 452 constexpr auto function = "boost::math::kurtosis_excess(const inverse_gamma_distribution<%1%>&)";453 454 RealType shape = dist.shape();455 RealType scale = dist.scale();456 457 RealType result = 0;458 if(false == detail::check_scale(function, scale, &result, Policy()))459 {460 return result;461 }462 if((shape <= 4) || !(boost::math::isfinite)(shape))463 {464 result = policies::raise_domain_error<RealType>(465 function,466 "Shape parameter is %1%, but for a defined kurtosis excess it must be > 4", shape, Policy());467 return result;468 }469 result = (30 * shape - 66) / ((shape - 3) * (shape - 4));470 return result;471}472 473template <class RealType, class Policy>474BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const inverse_gamma_distribution<RealType, Policy>& dist)475{476 constexpr auto function = "boost::math::kurtosis(const inverse_gamma_distribution<%1%>&)";477 RealType shape = dist.shape();478 RealType scale = dist.scale();479 480 RealType result = 0;481 482 if(false == detail::check_scale(function, scale, &result, Policy()))483 {484 return result;485 }486 if((shape <= 4) || !(boost::math::isfinite)(shape))487 {488 result = policies::raise_domain_error<RealType>(489 function,490 "Shape parameter is %1%, but for a defined kurtosis it must be > 4", shape, Policy());491 return result;492 }493 return kurtosis_excess(dist) + 3;494}495 496} // namespace math497} // namespace boost498 499// This include must be at the end, *after* the accessors500// for this distribution have been defined, in order to501// keep compilers that support two-phase lookup happy.502#include <boost/math/distributions/detail/derived_accessors.hpp>503 504#endif // BOOST_STATS_INVERSE_GAMMA_HPP505