ceph/src/boost/boost/math/distributions/weibull.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_WEIBULL_HPP
   7 #define BOOST_STATS_WEIBULL_HPP
   8
   9 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3668.htm
  10 // http://mathworld.wolfram.com/WeibullDistribution.html
  11
  12 #include <boost/math/distributions/fwd.hpp>
  13 #include <boost/math/special_functions/gamma.hpp>
  14 #include <boost/math/special_functions/log1p.hpp>
  15 #include <boost/math/special_functions/expm1.hpp>
  16 #include <boost/math/distributions/detail/common_error_handling.hpp>
  17 #include <boost/math/distributions/complement.hpp>
  18
  19 #include <utility>
  20
  21 namespace boost{ namespace math
  22 {
  23 namespace detail{
  24
  25 template <class RealType, class Policy>
  26 inline bool check_weibull_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_weibull_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_weibull(
  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_weibull_shape(function, shape, result, pol);
  65 }
  66
  67 } // namespace detail
  68
  69 template <class RealType = double, class Policy = policies::policy<> >
  70 class weibull_distribution
  71 {
  72 public:
  73    typedef RealType value_type;
  74    typedef Policy policy_type;
  75
  76    weibull_distribution(RealType l_shape, RealType l_scale = 1)
  77       : m_shape(l_shape), m_scale(l_scale)
  78    {
  79       RealType result;
  80       detail::check_weibull("boost::math::weibull_distribution<%1%>::weibull_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 typedef weibull_distribution<double> weibull;
 101
 102 template <class RealType, class Policy>
 103 inline const std::pair<RealType, RealType> range(const weibull_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 weibull_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    // A discontinuity at x == 0, so only support down to min_value.
 117 }
 118
 119 template <class RealType, class Policy>
 120 inline RealType pdf(const weibull_distribution<RealType, Policy>& dist, const RealType& x)
 121 {
 122    BOOST_MATH_STD_USING  // for ADL of std functions
 123
 124    static const char* function = "boost::math::pdf(const weibull_distribution<%1%>, %1%)";
 125
 126    RealType shape = dist.shape();
 127    RealType scale = dist.scale();
 128
 129    RealType result = 0;
 130    if(false == detail::check_weibull(function, scale, shape, &result, Policy()))
 131       return result;
 132    if(false == detail::check_weibull_x(function, x, &result, Policy()))
 133       return result;
 134
 135    if(x == 0)
 136    {
 137       if(shape == 1)
 138       {
 139          return 1 / scale;
 140       }
 141       if(shape > 1)
 142       {
 143          return 0;
 144       }
 145       return policies::raise_overflow_error<RealType>(function, 0, Policy());
 146    }
 147    result = exp(-pow(x / scale, shape));
 148    result *= pow(x / scale, shape - 1) * shape / scale;
 149
 150    return result;
 151 }
 152
 153 template <class RealType, class Policy>
 154 inline RealType cdf(const weibull_distribution<RealType, Policy>& dist, const RealType& x)
 155 {
 156    BOOST_MATH_STD_USING  // for ADL of std functions
 157
 158    static const char* function = "boost::math::cdf(const weibull_distribution<%1%>, %1%)";
 159
 160    RealType shape = dist.shape();
 161    RealType scale = dist.scale();
 162
 163    RealType result = 0;
 164    if(false == detail::check_weibull(function, scale, shape, &result, Policy()))
 165       return result;
 166    if(false == detail::check_weibull_x(function, x, &result, Policy()))
 167       return result;
 168
 169    result = -boost::math::expm1(-pow(x / scale, shape), Policy());
 170
 171    return result;
 172 }
 173
 174 template <class RealType, class Policy>
 175 inline RealType quantile(const weibull_distribution<RealType, Policy>& dist, const RealType& p)
 176 {
 177    BOOST_MATH_STD_USING  // for ADL of std functions
 178
 179    static const char* function = "boost::math::quantile(const weibull_distribution<%1%>, %1%)";
 180
 181    RealType shape = dist.shape();
 182    RealType scale = dist.scale();
 183
 184    RealType result = 0;
 185    if(false == detail::check_weibull(function, scale, shape, &result, Policy()))
 186       return result;
 187    if(false == detail::check_probability(function, p, &result, Policy()))
 188       return result;
 189
 190    if(p == 1)
 191       return policies::raise_overflow_error<RealType>(function, 0, Policy());
 192
 193    result = scale * pow(-boost::math::log1p(-p, Policy()), 1 / shape);
 194
 195    return result;
 196 }
 197
 198 template <class RealType, class Policy>
 199 inline RealType cdf(const complemented2_type<weibull_distribution<RealType, Policy>, RealType>& c)
 200 {
 201    BOOST_MATH_STD_USING  // for ADL of std functions
 202
 203    static const char* function = "boost::math::cdf(const weibull_distribution<%1%>, %1%)";
 204
 205    RealType shape = c.dist.shape();
 206    RealType scale = c.dist.scale();
 207
 208    RealType result = 0;
 209    if(false == detail::check_weibull(function, scale, shape, &result, Policy()))
 210       return result;
 211    if(false == detail::check_weibull_x(function, c.param, &result, Policy()))
 212       return result;
 213
 214    result = exp(-pow(c.param / scale, shape));
 215
 216    return result;
 217 }
 218
 219 template <class RealType, class Policy>
 220 inline RealType quantile(const complemented2_type<weibull_distribution<RealType, Policy>, RealType>& c)
 221 {
 222    BOOST_MATH_STD_USING  // for ADL of std functions
 223
 224    static const char* function = "boost::math::quantile(const weibull_distribution<%1%>, %1%)";
 225
 226    RealType shape = c.dist.shape();
 227    RealType scale = c.dist.scale();
 228    RealType q = c.param;
 229
 230    RealType result = 0;
 231    if(false == detail::check_weibull(function, scale, shape, &result, Policy()))
 232       return result;
 233    if(false == detail::check_probability(function, q, &result, Policy()))
 234       return result;
 235
 236    if(q == 0)
 237       return policies::raise_overflow_error<RealType>(function, 0, Policy());
 238
 239    result = scale * pow(-log(q), 1 / shape);
 240
 241    return result;
 242 }
 243
 244 template <class RealType, class Policy>
 245 inline RealType mean(const weibull_distribution<RealType, Policy>& dist)
 246 {
 247    BOOST_MATH_STD_USING  // for ADL of std functions
 248
 249    static const char* function = "boost::math::mean(const weibull_distribution<%1%>)";
 250
 251    RealType shape = dist.shape();
 252    RealType scale = dist.scale();
 253
 254    RealType result = 0;
 255    if(false == detail::check_weibull(function, scale, shape, &result, Policy()))
 256       return result;
 257
 258    result = scale * boost::math::tgamma(1 + 1 / shape, Policy());
 259    return result;
 260 }
 261
 262 template <class RealType, class Policy>
 263 inline RealType variance(const weibull_distribution<RealType, Policy>& dist)
 264 {
 265    RealType shape = dist.shape();
 266    RealType scale = dist.scale();
 267
 268    static const char* function = "boost::math::variance(const weibull_distribution<%1%>)";
 269
 270    RealType result = 0;
 271    if(false == detail::check_weibull(function, scale, shape, &result, Policy()))
 272    {
 273       return result;
 274    }
 275    result = boost::math::tgamma(1 + 1 / shape, Policy());
 276    result *= -result;
 277    result += boost::math::tgamma(1 + 2 / shape, Policy());
 278    result *= scale * scale;
 279    return result;
 280 }
 281
 282 template <class RealType, class Policy>
 283 inline RealType mode(const weibull_distribution<RealType, Policy>& dist)
 284 {
 285    BOOST_MATH_STD_USING  // for ADL of std function pow.
 286
 287    static const char* function = "boost::math::mode(const weibull_distribution<%1%>)";
 288
 289    RealType shape = dist.shape();
 290    RealType scale = dist.scale();
 291
 292    RealType result = 0;
 293    if(false == detail::check_weibull(function, scale, shape, &result, Policy()))
 294    {
 295       return result;
 296    }
 297    if(shape <= 1)
 298       return 0;
 299    result = scale * pow((shape - 1) / shape, 1 / shape);
 300    return result;
 301 }
 302
 303 template <class RealType, class Policy>
 304 inline RealType median(const weibull_distribution<RealType, Policy>& dist)
 305 {
 306    BOOST_MATH_STD_USING  // for ADL of std function pow.
 307
 308    static const char* function = "boost::math::median(const weibull_distribution<%1%>)";
 309
 310    RealType shape = dist.shape(); // Wikipedia k
 311    RealType scale = dist.scale(); // Wikipedia lambda
 312
 313    RealType result = 0;
 314    if(false == detail::check_weibull(function, scale, shape, &result, Policy()))
 315    {
 316       return result;
 317    }
 318    using boost::math::constants::ln_two;
 319    result = scale * pow(ln_two<RealType>(), 1 / shape);
 320    return result;
 321 }
 322
 323 template <class RealType, class Policy>
 324 inline RealType skewness(const weibull_distribution<RealType, Policy>& dist)
 325 {
 326    BOOST_MATH_STD_USING  // for ADL of std functions
 327
 328    static const char* function = "boost::math::skewness(const weibull_distribution<%1%>)";
 329
 330    RealType shape = dist.shape();
 331    RealType scale = dist.scale();
 332
 333    RealType result = 0;
 334    if(false == detail::check_weibull(function, scale, shape, &result, Policy()))
 335    {
 336       return result;
 337    }
 338    RealType g1, g2, g3, d;
 339
 340    g1 = boost::math::tgamma(1 + 1 / shape, Policy());
 341    g2 = boost::math::tgamma(1 + 2 / shape, Policy());
 342    g3 = boost::math::tgamma(1 + 3 / shape, Policy());
 343    d = pow(g2 - g1 * g1, RealType(1.5));
 344
 345    result = (2 * g1 * g1 * g1 - 3 * g1 * g2 + g3) / d;
 346    return result;
 347 }
 348
 349 template <class RealType, class Policy>
 350 inline RealType kurtosis_excess(const weibull_distribution<RealType, Policy>& dist)
 351 {
 352    BOOST_MATH_STD_USING  // for ADL of std functions
 353
 354    static const char* function = "boost::math::kurtosis_excess(const weibull_distribution<%1%>)";
 355
 356    RealType shape = dist.shape();
 357    RealType scale = dist.scale();
 358
 359    RealType result = 0;
 360    if(false == detail::check_weibull(function, scale, shape, &result, Policy()))
 361       return result;
 362
 363    RealType g1, g2, g3, g4, d, g1_2, g1_4;
 364
 365    g1 = boost::math::tgamma(1 + 1 / shape, Policy());
 366    g2 = boost::math::tgamma(1 + 2 / shape, Policy());
 367    g3 = boost::math::tgamma(1 + 3 / shape, Policy());
 368    g4 = boost::math::tgamma(1 + 4 / shape, Policy());
 369    g1_2 = g1 * g1;
 370    g1_4 = g1_2 * g1_2;
 371    d = g2 - g1_2;
 372    d *= d;
 373
 374    result = -6 * g1_4 + 12 * g1_2 * g2 - 3 * g2 * g2 - 4 * g1 * g3 + g4;
 375    result /= d;
 376    return result;
 377 }
 378
 379 template <class RealType, class Policy>
 380 inline RealType kurtosis(const weibull_distribution<RealType, Policy>& dist)
 381 {
 382    return kurtosis_excess(dist) + 3;
 383 }
 384
 385 } // namespace math
 386 } // namespace boost
 387
 388 // This include must be at the end, *after* the accessors
 389 // for this distribution have been defined, in order to
 390 // keep compilers that support two-phase lookup happy.
 391 #include <boost/math/distributions/detail/derived_accessors.hpp>
 392
 393 #endif // BOOST_STATS_WEIBULL_HPP
 394
 395