ceph/src/boost/boost/math/distributions/extreme_value.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_EXTREME_VALUE_HPP
   7 #define BOOST_STATS_EXTREME_VALUE_HPP
   8
   9 #include <boost/math/distributions/fwd.hpp>
  10 #include <boost/math/constants/constants.hpp>
  11 #include <boost/math/special_functions/log1p.hpp>
  12 #include <boost/math/special_functions/expm1.hpp>
  13 #include <boost/math/distributions/complement.hpp>
  14 #include <boost/math/distributions/detail/common_error_handling.hpp>
  15
  16 //
  17 // This is the maximum extreme value distribution, see
  18 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366g.htm
  19 // and http://mathworld.wolfram.com/ExtremeValueDistribution.html
  20 // Also known as a Fisher-Tippett distribution, a log-Weibull
  21 // distribution or a Gumbel distribution.
  22
  23 #include <utility>
  24 #include <cmath>
  25
  26 #ifdef _MSC_VER
  27 # pragma warning(push)
  28 # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
  29 #endif
  30
  31 namespace boost{ namespace math{
  32
  33 namespace detail{
  34 //
  35 // Error check:
  36 //
  37 template <class RealType, class Policy>
  38 inline bool verify_scale_b(const char* function, RealType b, RealType* presult, const Policy& pol)
  39 {
  40    if((b <= 0) || !(boost::math::isfinite)(b))
  41    {
  42       *presult = policies::raise_domain_error<RealType>(
  43          function,
  44          "The scale parameter \"b\" must be finite and > 0, but was: %1%.", b, pol);
  45       return false;
  46    }
  47    return true;
  48 }
  49
  50 } // namespace detail
  51
  52 template <class RealType = double, class Policy = policies::policy<> >
  53 class extreme_value_distribution
  54 {
  55 public:
  56    typedef RealType value_type;
  57    typedef Policy policy_type;
  58
  59    extreme_value_distribution(RealType a = 0, RealType b = 1)
  60       : m_a(a), m_b(b)
  61    {
  62       RealType err;
  63       detail::verify_scale_b("boost::math::extreme_value_distribution<%1%>::extreme_value_distribution", b, &err, Policy());
  64       detail::check_finite("boost::math::extreme_value_distribution<%1%>::extreme_value_distribution", a, &err, Policy());
  65    } // extreme_value_distribution
  66
  67    RealType location()const { return m_a; }
  68    RealType scale()const { return m_b; }
  69
  70 private:
  71    RealType m_a, m_b;
  72 };
  73
  74 typedef extreme_value_distribution<double> extreme_value;
  75
  76 #ifdef __cpp_deduction_guides
  77 template <class RealType>
  78 extreme_value_distribution(RealType)->extreme_value_distribution<typename boost::math::tools::promote_args<RealType>::type>;
  79 template <class RealType>
  80 extreme_value_distribution(RealType,RealType)->extreme_value_distribution<typename boost::math::tools::promote_args<RealType>::type>;
  81 #endif
  82
  83 template <class RealType, class Policy>
  84 inline const std::pair<RealType, RealType> range(const extreme_value_distribution<RealType, Policy>& /*dist*/)
  85 { // Range of permissible values for random variable x.
  86    using boost::math::tools::max_value;
  87    return std::pair<RealType, RealType>(
  88       std::numeric_limits<RealType>::has_infinity ? -std::numeric_limits<RealType>::infinity() : -max_value<RealType>(),
  89       std::numeric_limits<RealType>::has_infinity ? std::numeric_limits<RealType>::infinity() : max_value<RealType>());
  90 }
  91
  92 template <class RealType, class Policy>
  93 inline const std::pair<RealType, RealType> support(const extreme_value_distribution<RealType, Policy>& /*dist*/)
  94 { // Range of supported values for random variable x.
  95    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
  96    using boost::math::tools::max_value;
  97    return std::pair<RealType, RealType>(-max_value<RealType>(),  max_value<RealType>());
  98 }
  99
 100 template <class RealType, class Policy>
 101 inline RealType pdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x)
 102 {
 103    BOOST_MATH_STD_USING // for ADL of std functions
 104
 105    static const char* function = "boost::math::pdf(const extreme_value_distribution<%1%>&, %1%)";
 106
 107    RealType a = dist.location();
 108    RealType b = dist.scale();
 109    RealType result = 0;
 110    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
 111       return result;
 112    if(0 == detail::check_finite(function, a, &result, Policy()))
 113       return result;
 114    if((boost::math::isinf)(x))
 115       return 0.0f;
 116    if(0 == detail::check_x(function, x, &result, Policy()))
 117       return result;
 118    RealType e = (a - x) / b;
 119    if(e < tools::log_max_value<RealType>())
 120       result = exp(e) * exp(-exp(e)) / b;
 121    // else.... result *must* be zero since exp(e) is infinite...
 122    return result;
 123 } // pdf
 124
 125 template <class RealType, class Policy>
 126 inline RealType cdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x)
 127 {
 128    BOOST_MATH_STD_USING // for ADL of std functions
 129
 130    static const char* function = "boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)";
 131
 132    if((boost::math::isinf)(x))
 133       return x < 0 ? 0.0f : 1.0f;
 134    RealType a = dist.location();
 135    RealType b = dist.scale();
 136    RealType result = 0;
 137    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
 138       return result;
 139    if(0 == detail::check_finite(function, a, &result, Policy()))
 140       return result;
 141    if(0 == detail::check_finite(function, a, &result, Policy()))
 142       return result;
 143    if(0 == detail::check_x("boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)", x, &result, Policy()))
 144       return result;
 145
 146    result = exp(-exp((a-x)/b));
 147
 148    return result;
 149 } // cdf
 150
 151 template <class RealType, class Policy>
 152 RealType quantile(const extreme_value_distribution<RealType, Policy>& dist, const RealType& p)
 153 {
 154    BOOST_MATH_STD_USING // for ADL of std functions
 155
 156    static const char* function = "boost::math::quantile(const extreme_value_distribution<%1%>&, %1%)";
 157
 158    RealType a = dist.location();
 159    RealType b = dist.scale();
 160    RealType result = 0;
 161    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
 162       return result;
 163    if(0 == detail::check_finite(function, a, &result, Policy()))
 164       return result;
 165    if(0 == detail::check_probability(function, p, &result, Policy()))
 166       return result;
 167
 168    if(p == 0)
 169       return -policies::raise_overflow_error<RealType>(function, 0, Policy());
 170    if(p == 1)
 171       return policies::raise_overflow_error<RealType>(function, 0, Policy());
 172
 173    result = a - log(-log(p)) * b;
 174
 175    return result;
 176 } // quantile
 177
 178 template <class RealType, class Policy>
 179 inline RealType cdf(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c)
 180 {
 181    BOOST_MATH_STD_USING // for ADL of std functions
 182
 183    static const char* function = "boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)";
 184
 185    if((boost::math::isinf)(c.param))
 186       return c.param < 0 ? 1.0f : 0.0f;
 187    RealType a = c.dist.location();
 188    RealType b = c.dist.scale();
 189    RealType result = 0;
 190    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
 191       return result;
 192    if(0 == detail::check_finite(function, a, &result, Policy()))
 193       return result;
 194    if(0 == detail::check_x(function, c.param, &result, Policy()))
 195       return result;
 196
 197    result = -boost::math::expm1(-exp((a-c.param)/b), Policy());
 198
 199    return result;
 200 }
 201
 202 template <class RealType, class Policy>
 203 RealType quantile(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c)
 204 {
 205    BOOST_MATH_STD_USING // for ADL of std functions
 206
 207    static const char* function = "boost::math::quantile(const extreme_value_distribution<%1%>&, %1%)";
 208
 209    RealType a = c.dist.location();
 210    RealType b = c.dist.scale();
 211    RealType q = c.param;
 212    RealType result = 0;
 213    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
 214       return result;
 215    if(0 == detail::check_finite(function, a, &result, Policy()))
 216       return result;
 217    if(0 == detail::check_probability(function, q, &result, Policy()))
 218       return result;
 219
 220    if(q == 0)
 221       return policies::raise_overflow_error<RealType>(function, 0, Policy());
 222    if(q == 1)
 223       return -policies::raise_overflow_error<RealType>(function, 0, Policy());
 224
 225    result = a - log(-boost::math::log1p(-q, Policy())) * b;
 226
 227    return result;
 228 }
 229
 230 template <class RealType, class Policy>
 231 inline RealType mean(const extreme_value_distribution<RealType, Policy>& dist)
 232 {
 233    RealType a = dist.location();
 234    RealType b = dist.scale();
 235    RealType result = 0;
 236    if(0 == detail::verify_scale_b("boost::math::mean(const extreme_value_distribution<%1%>&)", b, &result, Policy()))
 237       return result;
 238    if (0 == detail::check_finite("boost::math::mean(const extreme_value_distribution<%1%>&)", a, &result, Policy()))
 239       return result;
 240    return a + constants::euler<RealType>() * b;
 241 }
 242
 243 template <class RealType, class Policy>
 244 inline RealType standard_deviation(const extreme_value_distribution<RealType, Policy>& dist)
 245 {
 246    BOOST_MATH_STD_USING // for ADL of std functions.
 247
 248    RealType b = dist.scale();
 249    RealType result = 0;
 250    if(0 == detail::verify_scale_b("boost::math::standard_deviation(const extreme_value_distribution<%1%>&)", b, &result, Policy()))
 251       return result;
 252    if(0 == detail::check_finite("boost::math::standard_deviation(const extreme_value_distribution<%1%>&)", dist.location(), &result, Policy()))
 253       return result;
 254    return constants::pi<RealType>() * b / sqrt(static_cast<RealType>(6));
 255 }
 256
 257 template <class RealType, class Policy>
 258 inline RealType mode(const extreme_value_distribution<RealType, Policy>& dist)
 259 {
 260    return dist.location();
 261 }
 262
 263 template <class RealType, class Policy>
 264 inline RealType median(const extreme_value_distribution<RealType, Policy>& dist)
 265 {
 266   using constants::ln_ln_two;
 267    return dist.location() - dist.scale() * ln_ln_two<RealType>();
 268 }
 269
 270 template <class RealType, class Policy>
 271 inline RealType skewness(const extreme_value_distribution<RealType, Policy>& /*dist*/)
 272 {
 273    //
 274    // This is 12 * sqrt(6) * zeta(3) / pi^3:
 275    // See http://mathworld.wolfram.com/ExtremeValueDistribution.html
 276    //
 277    return static_cast<RealType>(1.1395470994046486574927930193898461120875997958366L);
 278 }
 279
 280 template <class RealType, class Policy>
 281 inline RealType kurtosis(const extreme_value_distribution<RealType, Policy>& /*dist*/)
 282 {
 283    // See http://mathworld.wolfram.com/ExtremeValueDistribution.html
 284    return RealType(27) / 5;
 285 }
 286
 287 template <class RealType, class Policy>
 288 inline RealType kurtosis_excess(const extreme_value_distribution<RealType, Policy>& /*dist*/)
 289 {
 290    // See http://mathworld.wolfram.com/ExtremeValueDistribution.html
 291    return RealType(12) / 5;
 292 }
 293
 294
 295 } // namespace math
 296 } // namespace boost
 297
 298 #ifdef _MSC_VER
 299 # pragma warning(pop)
 300 #endif
 301
 302 // This include must be at the end, *after* the accessors
 303 // for this distribution have been defined, in order to
 304 // keep compilers that support two-phase lookup happy.
 305 #include <boost/math/distributions/detail/derived_accessors.hpp>
 306
 307 #endif // BOOST_STATS_EXTREME_VALUE_HPP