math/distributions/gamma.hpp

0001 //  Copyright John Maddock 2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_STATS_GAMMA_HPP
0007 #define BOOST_STATS_GAMMA_HPP
0008
0009 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
0010 // http://mathworld.wolfram.com/GammaDistribution.html
0011 // http://en.wikipedia.org/wiki/Gamma_distribution
0012
0013 #include <boost/math/distributions/fwd.hpp>
0014 #include <boost/math/special_functions/gamma.hpp>
0015 #include <boost/math/special_functions/digamma.hpp>
0016 #include <boost/math/distributions/detail/common_error_handling.hpp>
0017 #include <boost/math/distributions/complement.hpp>
0018
0019 #include <utility>
0020 #include <type_traits>
0021
0022 namespace boost{ namespace math
0023 {
0024 namespace detail
0025 {
0026
0027 template <class RealType, class Policy>
0028 inline bool check_gamma_shape(
0029       const char* function,
0030       RealType shape,
0031       RealType* result, const Policy& pol)
0032 {
0033    if((shape <= 0) || !(boost::math::isfinite)(shape))
0034    {
0035       *result = policies::raise_domain_error<RealType>(
0036          function,
0037          "Shape parameter is %1%, but must be > 0 !", shape, pol);
0038       return false;
0039    }
0040    return true;
0041 }
0042
0043 template <class RealType, class Policy>
0044 inline bool check_gamma_x(
0045       const char* function,
0046       RealType const& x,
0047       RealType* result, const Policy& pol)
0048 {
0049    if((x < 0) || !(boost::math::isfinite)(x))
0050    {
0051       *result = policies::raise_domain_error<RealType>(
0052          function,
0053          "Random variate is %1% but must be >= 0 !", x, pol);
0054       return false;
0055    }
0056    return true;
0057 }
0058
0059 template <class RealType, class Policy>
0060 inline bool check_gamma(
0061       const char* function,
0062       RealType scale,
0063       RealType shape,
0064       RealType* result, const Policy& pol)
0065 {
0066    return check_scale(function, scale, result, pol) && check_gamma_shape(function, shape, result, pol);
0067 }
0068
0069 } // namespace detail
0070
0071 template <class RealType = double, class Policy = policies::policy<> >
0072 class gamma_distribution
0073 {
0074 public:
0075    using value_type = RealType;
0076    using policy_type = Policy;
0077
0078    explicit gamma_distribution(RealType l_shape, RealType l_scale = 1)
0079       : m_shape(l_shape), m_scale(l_scale)
0080    {
0081       RealType result;
0082       detail::check_gamma("boost::math::gamma_distribution<%1%>::gamma_distribution", l_scale, l_shape, &result, Policy());
0083    }
0084
0085    RealType shape()const
0086    {
0087       return m_shape;
0088    }
0089
0090    RealType scale()const
0091    {
0092       return m_scale;
0093    }
0094 private:
0095    //
0096    // Data members:
0097    //
0098    RealType m_shape;     // distribution shape
0099    RealType m_scale;     // distribution scale
0100 };
0101
0102 // NO typedef because of clash with name of gamma function.
0103
0104 #ifdef __cpp_deduction_guides
0105 template <class RealType>
0106 gamma_distribution(RealType)->gamma_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0107 template <class RealType>
0108 gamma_distribution(RealType,RealType)->gamma_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0109 #endif
0110
0111 template <class RealType, class Policy>
0112 inline std::pair<RealType, RealType> range(const gamma_distribution<RealType, Policy>& /* dist */)
0113 { // Range of permissible values for random variable x.
0114    using boost::math::tools::max_value;
0115    return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
0116 }
0117
0118 template <class RealType, class Policy>
0119 inline std::pair<RealType, RealType> support(const gamma_distribution<RealType, Policy>& /* dist */)
0120 { // Range of supported values for random variable x.
0121    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0122    using boost::math::tools::max_value;
0123    using boost::math::tools::min_value;
0124    return std::pair<RealType, RealType>(min_value<RealType>(),  max_value<RealType>());
0125 }
0126
0127 template <class RealType, class Policy>
0128 inline RealType pdf(const gamma_distribution<RealType, Policy>& dist, const RealType& x)
0129 {
0130    BOOST_MATH_STD_USING  // for ADL of std functions
0131
0132    static const char* function = "boost::math::pdf(const gamma_distribution<%1%>&, %1%)";
0133
0134    RealType shape = dist.shape();
0135    RealType scale = dist.scale();
0136
0137    RealType result = 0;
0138    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
0139       return result;
0140    if(false == detail::check_gamma_x(function, x, &result, Policy()))
0141       return result;
0142
0143    if(x == 0)
0144    {
0145       return 0;
0146    }
0147    result = gamma_p_derivative(shape, x / scale, Policy()) / scale;
0148    return result;
0149 } // pdf
0150
0151 template <class RealType, class Policy>
0152 inline RealType logpdf(const gamma_distribution<RealType, Policy>& dist, const RealType& x)
0153 {
0154    BOOST_MATH_STD_USING  // for ADL of std functions
0155    using boost::math::lgamma;
0156
0157    static const char* function = "boost::math::logpdf(const gamma_distribution<%1%>&, %1%)";
0158
0159    RealType k = dist.shape();
0160    RealType theta = dist.scale();
0161
0162    RealType result = -std::numeric_limits<RealType>::infinity();
0163    if(false == detail::check_gamma(function, theta, k, &result, Policy()))
0164       return result;
0165    if(false == detail::check_gamma_x(function, x, &result, Policy()))
0166       return result;
0167
0168    if(x == 0)
0169    {
0170       return std::numeric_limits<RealType>::quiet_NaN();
0171    }
0172
0173    result = -k*log(theta) + (k-1)*log(x) - lgamma(k) - (x/theta);
0174
0175    return result;
0176 } // logpdf
0177
0178 template <class RealType, class Policy>
0179 inline RealType cdf(const gamma_distribution<RealType, Policy>& dist, const RealType& x)
0180 {
0181    BOOST_MATH_STD_USING  // for ADL of std functions
0182
0183    static const char* function = "boost::math::cdf(const gamma_distribution<%1%>&, %1%)";
0184
0185    RealType shape = dist.shape();
0186    RealType scale = dist.scale();
0187
0188    RealType result = 0;
0189    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
0190       return result;
0191    if(false == detail::check_gamma_x(function, x, &result, Policy()))
0192       return result;
0193
0194    result = boost::math::gamma_p(shape, x / scale, Policy());
0195    return result;
0196 } // cdf
0197
0198 template <class RealType, class Policy>
0199 inline RealType quantile(const gamma_distribution<RealType, Policy>& dist, const RealType& p)
0200 {
0201    BOOST_MATH_STD_USING  // for ADL of std functions
0202
0203    static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";
0204
0205    RealType shape = dist.shape();
0206    RealType scale = dist.scale();
0207
0208    RealType result = 0;
0209    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
0210       return result;
0211    if(false == detail::check_probability(function, p, &result, Policy()))
0212       return result;
0213
0214    if(p == 1)
0215       return policies::raise_overflow_error<RealType>(function, 0, Policy());
0216
0217    result = gamma_p_inv(shape, p, Policy()) * scale;
0218
0219    return result;
0220 }
0221
0222 template <class RealType, class Policy>
0223 inline RealType cdf(const complemented2_type<gamma_distribution<RealType, Policy>, RealType>& c)
0224 {
0225    BOOST_MATH_STD_USING  // for ADL of std functions
0226
0227    static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";
0228
0229    RealType shape = c.dist.shape();
0230    RealType scale = c.dist.scale();
0231
0232    RealType result = 0;
0233    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
0234       return result;
0235    if(false == detail::check_gamma_x(function, c.param, &result, Policy()))
0236       return result;
0237
0238    result = gamma_q(shape, c.param / scale, Policy());
0239
0240    return result;
0241 }
0242
0243 template <class RealType, class Policy>
0244 inline RealType quantile(const complemented2_type<gamma_distribution<RealType, Policy>, RealType>& c)
0245 {
0246    BOOST_MATH_STD_USING  // for ADL of std functions
0247
0248    static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";
0249
0250    RealType shape = c.dist.shape();
0251    RealType scale = c.dist.scale();
0252    RealType q = c.param;
0253
0254    RealType result = 0;
0255    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
0256       return result;
0257    if(false == detail::check_probability(function, q, &result, Policy()))
0258       return result;
0259
0260    if(q == 0)
0261       return policies::raise_overflow_error<RealType>(function, 0, Policy());
0262
0263    result = gamma_q_inv(shape, q, Policy()) * scale;
0264
0265    return result;
0266 }
0267
0268 template <class RealType, class Policy>
0269 inline RealType mean(const gamma_distribution<RealType, Policy>& dist)
0270 {
0271    BOOST_MATH_STD_USING  // for ADL of std functions
0272
0273    static const char* function = "boost::math::mean(const gamma_distribution<%1%>&)";
0274
0275    RealType shape = dist.shape();
0276    RealType scale = dist.scale();
0277
0278    RealType result = 0;
0279    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
0280       return result;
0281
0282    result = shape * scale;
0283    return result;
0284 }
0285
0286 template <class RealType, class Policy>
0287 inline RealType variance(const gamma_distribution<RealType, Policy>& dist)
0288 {
0289    BOOST_MATH_STD_USING  // for ADL of std functions
0290
0291    static const char* function = "boost::math::variance(const gamma_distribution<%1%>&)";
0292
0293    RealType shape = dist.shape();
0294    RealType scale = dist.scale();
0295
0296    RealType result = 0;
0297    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
0298       return result;
0299
0300    result = shape * scale * scale;
0301    return result;
0302 }
0303
0304 template <class RealType, class Policy>
0305 inline RealType mode(const gamma_distribution<RealType, Policy>& dist)
0306 {
0307    BOOST_MATH_STD_USING  // for ADL of std functions
0308
0309    static const char* function = "boost::math::mode(const gamma_distribution<%1%>&)";
0310
0311    RealType shape = dist.shape();
0312    RealType scale = dist.scale();
0313
0314    RealType result = 0;
0315    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
0316       return result;
0317
0318    if(shape < 1)
0319       return policies::raise_domain_error<RealType>(
0320          function,
0321          "The mode of the gamma distribution is only defined for values of the shape parameter >= 1, but got %1%.",
0322          shape, Policy());
0323
0324    result = (shape - 1) * scale;
0325    return result;
0326 }
0327
0328 //template <class RealType, class Policy>
0329 //inline RealType median(const gamma_distribution<RealType, Policy>& dist)
0330 //{  // Rely on default definition in derived accessors.
0331 //}
0332
0333 template <class RealType, class Policy>
0334 inline RealType skewness(const gamma_distribution<RealType, Policy>& dist)
0335 {
0336    BOOST_MATH_STD_USING  // for ADL of std functions
0337
0338    static const char* function = "boost::math::skewness(const gamma_distribution<%1%>&)";
0339
0340    RealType shape = dist.shape();
0341    RealType scale = dist.scale();
0342
0343    RealType result = 0;
0344    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
0345       return result;
0346
0347    result = 2 / sqrt(shape);
0348    return result;
0349 }
0350
0351 template <class RealType, class Policy>
0352 inline RealType kurtosis_excess(const gamma_distribution<RealType, Policy>& dist)
0353 {
0354    BOOST_MATH_STD_USING  // for ADL of std functions
0355
0356    static const char* function = "boost::math::kurtosis_excess(const gamma_distribution<%1%>&)";
0357
0358    RealType shape = dist.shape();
0359    RealType scale = dist.scale();
0360
0361    RealType result = 0;
0362    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
0363       return result;
0364
0365    result = 6 / shape;
0366    return result;
0367 }
0368
0369 template <class RealType, class Policy>
0370 inline RealType kurtosis(const gamma_distribution<RealType, Policy>& dist)
0371 {
0372    return kurtosis_excess(dist) + 3;
0373 }
0374
0375 template <class RealType, class Policy>
0376 inline RealType entropy(const gamma_distribution<RealType, Policy>& dist)
0377 {
0378    RealType k = dist.shape();
0379    RealType theta = dist.scale();
0380    using std::log;
0381    return k + log(theta) + lgamma(k) + (1-k)*digamma(k);
0382 }
0383
0384 } // namespace math
0385 } // namespace boost
0386
0387 // This include must be at the end, *after* the accessors
0388 // for this distribution have been defined, in order to
0389 // keep compilers that support two-phase lookup happy.
0390 #include <boost/math/distributions/detail/derived_accessors.hpp>
0391
0392 #endif // BOOST_STATS_GAMMA_HPP
0393
0394