ceph/src/boost/boost/math/distributions/gamma.hpp

   1 //  Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_STATS_GAMMA_HPP
   7 #define BOOST_STATS_GAMMA_HPP
   8
   9 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
  10 // http://mathworld.wolfram.com/GammaDistribution.html
  11 // http://en.wikipedia.org/wiki/Gamma_distribution
  12
  13 #include <boost/math/distributions/fwd.hpp>
  14 #include <boost/math/special_functions/gamma.hpp>
  15 #include <boost/math/distributions/detail/common_error_handling.hpp>
  16 #include <boost/math/distributions/complement.hpp>
  17
  18 #include <utility>
  19
  20 namespace boost{ namespace math
  21 {
  22 namespace detail
  23 {
  24
  25 template <class RealType, class Policy>
  26 inline bool check_gamma_shape(
  27       const char* function,
  28       RealType shape,
  29       RealType* result, const Policy& pol)
  30 {
  31    if((shape <= 0) || !(boost::math::isfinite)(shape))
  32    {
  33       *result = policies::raise_domain_error<RealType>(
  34          function,
  35          "Shape parameter is %1%, but must be > 0 !", shape, pol);
  36       return false;
  37    }
  38    return true;
  39 }
  40
  41 template <class RealType, class Policy>
  42 inline bool check_gamma_x(
  43       const char* function,
  44       RealType const& x,
  45       RealType* result, const Policy& pol)
  46 {
  47    if((x < 0) || !(boost::math::isfinite)(x))
  48    {
  49       *result = policies::raise_domain_error<RealType>(
  50          function,
  51          "Random variate is %1% but must be >= 0 !", x, pol);
  52       return false;
  53    }
  54    return true;
  55 }
  56
  57 template <class RealType, class Policy>
  58 inline bool check_gamma(
  59       const char* function,
  60       RealType scale,
  61       RealType shape,
  62       RealType* result, const Policy& pol)
  63 {
  64    return check_scale(function, scale, result, pol) && check_gamma_shape(function, shape, result, pol);
  65 }
  66
  67 } // namespace detail
  68
  69 template <class RealType = double, class Policy = policies::policy<> >
  70 class gamma_distribution
  71 {
  72 public:
  73    typedef RealType value_type;
  74    typedef Policy policy_type;
  75
  76    gamma_distribution(RealType l_shape, RealType l_scale = 1)
  77       : m_shape(l_shape), m_scale(l_scale)
  78    {
  79       RealType result;
  80       detail::check_gamma("boost::math::gamma_distribution<%1%>::gamma_distribution", l_scale, l_shape, &result, Policy());
  81    }
  82
  83    RealType shape()const
  84    {
  85       return m_shape;
  86    }
  87
  88    RealType scale()const
  89    {
  90       return m_scale;
  91    }
  92 private:
  93    //
  94    // Data members:
  95    //
  96    RealType m_shape;     // distribution shape
  97    RealType m_scale;     // distribution scale
  98 };
  99
 100 // NO typedef because of clash with name of gamma function.
 101
 102 template <class RealType, class Policy>
 103 inline const std::pair<RealType, RealType> range(const gamma_distribution<RealType, Policy>& /* dist */)
 104 { // Range of permissible values for random variable x.
 105    using boost::math::tools::max_value;
 106    return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
 107 }
 108
 109 template <class RealType, class Policy>
 110 inline const std::pair<RealType, RealType> support(const gamma_distribution<RealType, Policy>& /* dist */)
 111 { // Range of supported values for random variable x.
 112    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
 113    using boost::math::tools::max_value;
 114    using boost::math::tools::min_value;
 115    return std::pair<RealType, RealType>(min_value<RealType>(),  max_value<RealType>());
 116 }
 117
 118 template <class RealType, class Policy>
 119 inline RealType pdf(const gamma_distribution<RealType, Policy>& dist, const RealType& x)
 120 {
 121    BOOST_MATH_STD_USING  // for ADL of std functions
 122
 123    static const char* function = "boost::math::pdf(const gamma_distribution<%1%>&, %1%)";
 124
 125    RealType shape = dist.shape();
 126    RealType scale = dist.scale();
 127
 128    RealType result = 0;
 129    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
 130       return result;
 131    if(false == detail::check_gamma_x(function, x, &result, Policy()))
 132       return result;
 133
 134    if(x == 0)
 135    {
 136       return 0;
 137    }
 138    result = gamma_p_derivative(shape, x / scale, Policy()) / scale;
 139    return result;
 140 } // pdf
 141
 142 template <class RealType, class Policy>
 143 inline RealType cdf(const gamma_distribution<RealType, Policy>& dist, const RealType& x)
 144 {
 145    BOOST_MATH_STD_USING  // for ADL of std functions
 146
 147    static const char* function = "boost::math::cdf(const gamma_distribution<%1%>&, %1%)";
 148
 149    RealType shape = dist.shape();
 150    RealType scale = dist.scale();
 151
 152    RealType result = 0;
 153    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
 154       return result;
 155    if(false == detail::check_gamma_x(function, x, &result, Policy()))
 156       return result;
 157
 158    result = boost::math::gamma_p(shape, x / scale, Policy());
 159    return result;
 160 } // cdf
 161
 162 template <class RealType, class Policy>
 163 inline RealType quantile(const gamma_distribution<RealType, Policy>& dist, const RealType& p)
 164 {
 165    BOOST_MATH_STD_USING  // for ADL of std functions
 166
 167    static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";
 168
 169    RealType shape = dist.shape();
 170    RealType scale = dist.scale();
 171
 172    RealType result = 0;
 173    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
 174       return result;
 175    if(false == detail::check_probability(function, p, &result, Policy()))
 176       return result;
 177
 178    if(p == 1)
 179       return policies::raise_overflow_error<RealType>(function, 0, Policy());
 180
 181    result = gamma_p_inv(shape, p, Policy()) * scale;
 182
 183    return result;
 184 }
 185
 186 template <class RealType, class Policy>
 187 inline RealType cdf(const complemented2_type<gamma_distribution<RealType, Policy>, RealType>& c)
 188 {
 189    BOOST_MATH_STD_USING  // for ADL of std functions
 190
 191    static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";
 192
 193    RealType shape = c.dist.shape();
 194    RealType scale = c.dist.scale();
 195
 196    RealType result = 0;
 197    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
 198       return result;
 199    if(false == detail::check_gamma_x(function, c.param, &result, Policy()))
 200       return result;
 201
 202    result = gamma_q(shape, c.param / scale, Policy());
 203
 204    return result;
 205 }
 206
 207 template <class RealType, class Policy>
 208 inline RealType quantile(const complemented2_type<gamma_distribution<RealType, Policy>, RealType>& c)
 209 {
 210    BOOST_MATH_STD_USING  // for ADL of std functions
 211
 212    static const char* function = "boost::math::quantile(const gamma_distribution<%1%>&, %1%)";
 213
 214    RealType shape = c.dist.shape();
 215    RealType scale = c.dist.scale();
 216    RealType q = c.param;
 217
 218    RealType result = 0;
 219    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
 220       return result;
 221    if(false == detail::check_probability(function, q, &result, Policy()))
 222       return result;
 223
 224    if(q == 0)
 225       return policies::raise_overflow_error<RealType>(function, 0, Policy());
 226
 227    result = gamma_q_inv(shape, q, Policy()) * scale;
 228
 229    return result;
 230 }
 231
 232 template <class RealType, class Policy>
 233 inline RealType mean(const gamma_distribution<RealType, Policy>& dist)
 234 {
 235    BOOST_MATH_STD_USING  // for ADL of std functions
 236
 237    static const char* function = "boost::math::mean(const gamma_distribution<%1%>&)";
 238
 239    RealType shape = dist.shape();
 240    RealType scale = dist.scale();
 241
 242    RealType result = 0;
 243    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
 244       return result;
 245
 246    result = shape * scale;
 247    return result;
 248 }
 249
 250 template <class RealType, class Policy>
 251 inline RealType variance(const gamma_distribution<RealType, Policy>& dist)
 252 {
 253    BOOST_MATH_STD_USING  // for ADL of std functions
 254
 255    static const char* function = "boost::math::variance(const gamma_distribution<%1%>&)";
 256
 257    RealType shape = dist.shape();
 258    RealType scale = dist.scale();
 259
 260    RealType result = 0;
 261    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
 262       return result;
 263
 264    result = shape * scale * scale;
 265    return result;
 266 }
 267
 268 template <class RealType, class Policy>
 269 inline RealType mode(const gamma_distribution<RealType, Policy>& dist)
 270 {
 271    BOOST_MATH_STD_USING  // for ADL of std functions
 272
 273    static const char* function = "boost::math::mode(const gamma_distribution<%1%>&)";
 274
 275    RealType shape = dist.shape();
 276    RealType scale = dist.scale();
 277
 278    RealType result = 0;
 279    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
 280       return result;
 281
 282    if(shape < 1)
 283       return policies::raise_domain_error<RealType>(
 284          function,
 285          "The mode of the gamma distribution is only defined for values of the shape parameter >= 1, but got %1%.",
 286          shape, Policy());
 287
 288    result = (shape - 1) * scale;
 289    return result;
 290 }
 291
 292 //template <class RealType, class Policy>
 293 //inline RealType median(const gamma_distribution<RealType, Policy>& dist)
 294 //{  // Rely on default definition in derived accessors.
 295 //}
 296
 297 template <class RealType, class Policy>
 298 inline RealType skewness(const gamma_distribution<RealType, Policy>& dist)
 299 {
 300    BOOST_MATH_STD_USING  // for ADL of std functions
 301
 302    static const char* function = "boost::math::skewness(const gamma_distribution<%1%>&)";
 303
 304    RealType shape = dist.shape();
 305    RealType scale = dist.scale();
 306
 307    RealType result = 0;
 308    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
 309       return result;
 310
 311    result = 2 / sqrt(shape);
 312    return result;
 313 }
 314
 315 template <class RealType, class Policy>
 316 inline RealType kurtosis_excess(const gamma_distribution<RealType, Policy>& dist)
 317 {
 318    BOOST_MATH_STD_USING  // for ADL of std functions
 319
 320    static const char* function = "boost::math::kurtosis_excess(const gamma_distribution<%1%>&)";
 321
 322    RealType shape = dist.shape();
 323    RealType scale = dist.scale();
 324
 325    RealType result = 0;
 326    if(false == detail::check_gamma(function, scale, shape, &result, Policy()))
 327       return result;
 328
 329    result = 6 / shape;
 330    return result;
 331 }
 332
 333 template <class RealType, class Policy>
 334 inline RealType kurtosis(const gamma_distribution<RealType, Policy>& dist)
 335 {
 336    return kurtosis_excess(dist) + 3;
 337 }
 338
 339 } // namespace math
 340 } // namespace boost
 341
 342 // This include must be at the end, *after* the accessors
 343 // for this distribution have been defined, in order to
 344 // keep compilers that support two-phase lookup happy.
 345 #include <boost/math/distributions/detail/derived_accessors.hpp>
 346
 347 #endif // BOOST_STATS_GAMMA_HPP
 348
 349