397 lines · plain
1// Copyright John Maddock 2006.2// Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_STATS_GAMMA_HPP8#define BOOST_STATS_GAMMA_HPP9 10// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm11// http://mathworld.wolfram.com/GammaDistribution.html12// http://en.wikipedia.org/wiki/Gamma_distribution13 14#include <boost/math/tools/config.hpp>15#include <boost/math/tools/tuple.hpp>16#include <boost/math/tools/numeric_limits.hpp>17#include <boost/math/distributions/fwd.hpp>18#include <boost/math/special_functions/gamma.hpp>19#include <boost/math/special_functions/digamma.hpp>20#include <boost/math/distributions/detail/common_error_handling.hpp>21#include <boost/math/distributions/complement.hpp>22 23namespace boost{ namespace math24{25namespace detail26{27 28template <class RealType, class Policy>29BOOST_MATH_GPU_ENABLED inline bool check_gamma_shape(30 const char* function,31 RealType shape,32 RealType* result, const Policy& pol)33{34 if((shape <= 0) || !(boost::math::isfinite)(shape))35 {36 *result = policies::raise_domain_error<RealType>(37 function,38 "Shape parameter is %1%, but must be > 0 !", shape, pol);39 return false;40 }41 return true;42}43 44template <class RealType, class Policy>45BOOST_MATH_GPU_ENABLED inline bool check_gamma_x(46 const char* function,47 RealType const& x,48 RealType* result, const Policy& pol)49{50 if((x < 0) || !(boost::math::isfinite)(x))51 {52 *result = policies::raise_domain_error<RealType>(53 function,54 "Random variate is %1% but must be >= 0 !", x, pol);55 return false;56 }57 return true;58}59 60template <class RealType, class Policy>61BOOST_MATH_GPU_ENABLED inline bool check_gamma(62 const char* function,63 RealType scale,64 RealType shape,65 RealType* result, const Policy& pol)66{67 return check_scale(function, scale, result, pol) && check_gamma_shape(function, shape, result, pol);68}69 70} // namespace detail71 72template <class RealType = double, class Policy = policies::policy<> >73class gamma_distribution74{75public:76 using value_type = RealType;77 using policy_type = Policy;78 79 BOOST_MATH_GPU_ENABLED explicit gamma_distribution(RealType l_shape, RealType l_scale = 1)80 : m_shape(l_shape), m_scale(l_scale)81 {82 RealType result;83 detail::check_gamma("boost::math::gamma_distribution<%1%>::gamma_distribution", l_scale, l_shape, &result, Policy());84 }85 86 BOOST_MATH_GPU_ENABLED RealType shape()const87 {88 return m_shape;89 }90 91 BOOST_MATH_GPU_ENABLED RealType scale()const92 {93 return m_scale;94 }95private:96 //97 // Data members:98 //99 RealType m_shape; // distribution shape100 RealType m_scale; // distribution scale101};102 103// NO typedef because of clash with name of gamma function.104 105#ifdef __cpp_deduction_guides106template <class RealType>107gamma_distribution(RealType)->gamma_distribution<typename boost::math::tools::promote_args<RealType>::type>;108template <class RealType>109gamma_distribution(RealType,RealType)->gamma_distribution<typename boost::math::tools::promote_args<RealType>::type>;110#endif111 112template <class RealType, class Policy>113BOOST_MATH_GPU_ENABLED inline boost::math::pair<RealType, RealType> range(const gamma_distribution<RealType, Policy>& /* dist */)114{ // Range of permissible values for random variable x.115 using boost::math::tools::max_value;116 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());117}118 119template <class RealType, class Policy>120BOOST_MATH_GPU_ENABLED inline boost::math::pair<RealType, RealType> support(const gamma_distribution<RealType, Policy>& /* dist */)121{ // Range of supported values for random variable x.122 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.123 using boost::math::tools::max_value;124 using boost::math::tools::min_value;125 return boost::math::pair<RealType, RealType>(min_value<RealType>(), max_value<RealType>());126}127 128template <class RealType, class Policy>129BOOST_MATH_GPU_ENABLED inline RealType pdf(const gamma_distribution<RealType, Policy>& dist, const RealType& x)130{131 BOOST_MATH_STD_USING // for ADL of std functions132 133 constexpr auto function = "boost::math::pdf(const gamma_distribution<%1%>&, %1%)";134 135 RealType shape = dist.shape();136 RealType scale = dist.scale();137 138 RealType result = 0;139 if(false == detail::check_gamma(function, scale, shape, &result, Policy()))140 return result;141 if(false == detail::check_gamma_x(function, x, &result, Policy()))142 return result;143 144 if(x == 0)145 {146 return 0;147 }148 result = gamma_p_derivative(shape, x / scale, Policy()) / scale;149 return result;150} // pdf151 152template <class RealType, class Policy>153BOOST_MATH_GPU_ENABLED inline RealType logpdf(const gamma_distribution<RealType, Policy>& dist, const RealType& x)154{155 BOOST_MATH_STD_USING // for ADL of std functions156 using boost::math::lgamma;157 158 constexpr auto function = "boost::math::logpdf(const gamma_distribution<%1%>&, %1%)";159 160 RealType k = dist.shape();161 RealType theta = dist.scale();162 163 RealType result = -boost::math::numeric_limits<RealType>::infinity();164 if(false == detail::check_gamma(function, theta, k, &result, Policy()))165 return result;166 if(false == detail::check_gamma_x(function, x, &result, Policy()))167 return result;168 169 if(x == 0)170 {171 return boost::math::numeric_limits<RealType>::quiet_NaN();172 }173 174 result = -k*log(theta) + (k-1)*log(x) - lgamma(k) - (x/theta);175 176 return result;177} // logpdf178 179template <class RealType, class Policy>180BOOST_MATH_GPU_ENABLED inline RealType cdf(const gamma_distribution<RealType, Policy>& dist, const RealType& x)181{182 BOOST_MATH_STD_USING // for ADL of std functions183 184 constexpr auto function = "boost::math::cdf(const gamma_distribution<%1%>&, %1%)";185 186 RealType shape = dist.shape();187 RealType scale = dist.scale();188 189 RealType result = 0;190 if(false == detail::check_gamma(function, scale, shape, &result, Policy()))191 return result;192 if(false == detail::check_gamma_x(function, x, &result, Policy()))193 return result;194 195 result = boost::math::gamma_p(shape, x / scale, Policy());196 return result;197} // cdf198 199template <class RealType, class Policy>200BOOST_MATH_GPU_ENABLED inline RealType quantile(const gamma_distribution<RealType, Policy>& dist, const RealType& p)201{202 BOOST_MATH_STD_USING // for ADL of std functions203 204 constexpr auto function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";205 206 RealType shape = dist.shape();207 RealType scale = dist.scale();208 209 RealType result = 0;210 if(false == detail::check_gamma(function, scale, shape, &result, Policy()))211 return result;212 if(false == detail::check_probability(function, p, &result, Policy()))213 return result;214 215 if(p == 1)216 return policies::raise_overflow_error<RealType>(function, 0, Policy());217 218 result = gamma_p_inv(shape, p, Policy()) * scale;219 220 return result;221}222 223template <class RealType, class Policy>224BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<gamma_distribution<RealType, Policy>, RealType>& c)225{226 BOOST_MATH_STD_USING // for ADL of std functions227 228 constexpr auto function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";229 230 RealType shape = c.dist.shape();231 RealType scale = c.dist.scale();232 233 RealType result = 0;234 if(false == detail::check_gamma(function, scale, shape, &result, Policy()))235 return result;236 if(false == detail::check_gamma_x(function, c.param, &result, Policy()))237 return result;238 239 result = gamma_q(shape, c.param / scale, Policy());240 241 return result;242}243 244template <class RealType, class Policy>245BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<gamma_distribution<RealType, Policy>, RealType>& c)246{247 BOOST_MATH_STD_USING // for ADL of std functions248 249 constexpr auto function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";250 251 RealType shape = c.dist.shape();252 RealType scale = c.dist.scale();253 RealType q = c.param;254 255 RealType result = 0;256 if(false == detail::check_gamma(function, scale, shape, &result, Policy()))257 return result;258 if(false == detail::check_probability(function, q, &result, Policy()))259 return result;260 261 if(q == 0)262 return policies::raise_overflow_error<RealType>(function, 0, Policy());263 264 result = gamma_q_inv(shape, q, Policy()) * scale;265 266 return result;267}268 269template <class RealType, class Policy>270BOOST_MATH_GPU_ENABLED inline RealType mean(const gamma_distribution<RealType, Policy>& dist)271{272 BOOST_MATH_STD_USING // for ADL of std functions273 274 constexpr auto function = "boost::math::mean(const gamma_distribution<%1%>&)";275 276 RealType shape = dist.shape();277 RealType scale = dist.scale();278 279 RealType result = 0;280 if(false == detail::check_gamma(function, scale, shape, &result, Policy()))281 return result;282 283 result = shape * scale;284 return result;285}286 287template <class RealType, class Policy>288BOOST_MATH_GPU_ENABLED inline RealType variance(const gamma_distribution<RealType, Policy>& dist)289{290 BOOST_MATH_STD_USING // for ADL of std functions291 292 constexpr auto function = "boost::math::variance(const gamma_distribution<%1%>&)";293 294 RealType shape = dist.shape();295 RealType scale = dist.scale();296 297 RealType result = 0;298 if(false == detail::check_gamma(function, scale, shape, &result, Policy()))299 return result;300 301 result = shape * scale * scale;302 return result;303}304 305template <class RealType, class Policy>306BOOST_MATH_GPU_ENABLED inline RealType mode(const gamma_distribution<RealType, Policy>& dist)307{308 BOOST_MATH_STD_USING // for ADL of std functions309 310 constexpr auto function = "boost::math::mode(const gamma_distribution<%1%>&)";311 312 RealType shape = dist.shape();313 RealType scale = dist.scale();314 315 RealType result = 0;316 if(false == detail::check_gamma(function, scale, shape, &result, Policy()))317 return result;318 319 if(shape < 1)320 return policies::raise_domain_error<RealType>(321 function,322 "The mode of the gamma distribution is only defined for values of the shape parameter >= 1, but got %1%.",323 shape, Policy());324 325 result = (shape - 1) * scale;326 return result;327}328 329//template <class RealType, class Policy>330//inline RealType median(const gamma_distribution<RealType, Policy>& dist)331//{ // Rely on default definition in derived accessors.332//}333 334template <class RealType, class Policy>335BOOST_MATH_GPU_ENABLED inline RealType skewness(const gamma_distribution<RealType, Policy>& dist)336{337 BOOST_MATH_STD_USING // for ADL of std functions338 339 constexpr auto function = "boost::math::skewness(const gamma_distribution<%1%>&)";340 341 RealType shape = dist.shape();342 RealType scale = dist.scale();343 344 RealType result = 0;345 if(false == detail::check_gamma(function, scale, shape, &result, Policy()))346 return result;347 348 result = 2 / sqrt(shape);349 return result;350}351 352template <class RealType, class Policy>353BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const gamma_distribution<RealType, Policy>& dist)354{355 BOOST_MATH_STD_USING // for ADL of std functions356 357 constexpr auto function = "boost::math::kurtosis_excess(const gamma_distribution<%1%>&)";358 359 RealType shape = dist.shape();360 RealType scale = dist.scale();361 362 RealType result = 0;363 if(false == detail::check_gamma(function, scale, shape, &result, Policy()))364 return result;365 366 result = 6 / shape;367 return result;368}369 370template <class RealType, class Policy>371BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const gamma_distribution<RealType, Policy>& dist)372{373 return kurtosis_excess(dist) + 3;374}375 376template <class RealType, class Policy>377BOOST_MATH_GPU_ENABLED inline RealType entropy(const gamma_distribution<RealType, Policy>& dist)378{379 BOOST_MATH_STD_USING380 381 RealType k = dist.shape();382 RealType theta = dist.scale();383 return k + log(theta) + boost::math::lgamma(k) + (1-k)*digamma(k);384}385 386} // namespace math387} // namespace boost388 389// This include must be at the end, *after* the accessors390// for this distribution have been defined, in order to391// keep compilers that support two-phase lookup happy.392#include <boost/math/distributions/detail/derived_accessors.hpp>393 394#endif // BOOST_STATS_GAMMA_HPP395 396 397