math/distributions/extreme_value.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_EXTREME_VALUE_HPP
0007 #define BOOST_STATS_EXTREME_VALUE_HPP
0008
0009 #include <boost/math/distributions/fwd.hpp>
0010 #include <boost/math/constants/constants.hpp>
0011 #include <boost/math/special_functions/log1p.hpp>
0012 #include <boost/math/special_functions/expm1.hpp>
0013 #include <boost/math/distributions/complement.hpp>
0014 #include <boost/math/distributions/detail/common_error_handling.hpp>
0015
0016 //
0017 // This is the maximum extreme value distribution, see
0018 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366g.htm
0019 // and http://mathworld.wolfram.com/ExtremeValueDistribution.html
0020 // Also known as a Fisher-Tippett distribution, a log-Weibull
0021 // distribution or a Gumbel distribution.
0022
0023 #include <utility>
0024 #include <cmath>
0025
0026 #ifdef _MSC_VER
0027 # pragma warning(push)
0028 # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
0029 #endif
0030
0031 namespace boost{ namespace math{
0032
0033 namespace detail{
0034 //
0035 // Error check:
0036 //
0037 template <class RealType, class Policy>
0038 inline bool verify_scale_b(const char* function, RealType b, RealType* presult, const Policy& pol)
0039 {
0040    if((b <= 0) || !(boost::math::isfinite)(b))
0041    {
0042       *presult = policies::raise_domain_error<RealType>(
0043          function,
0044          "The scale parameter \"b\" must be finite and > 0, but was: %1%.", b, pol);
0045       return false;
0046    }
0047    return true;
0048 }
0049
0050 } // namespace detail
0051
0052 template <class RealType = double, class Policy = policies::policy<> >
0053 class extreme_value_distribution
0054 {
0055 public:
0056    using value_type = RealType;
0057    using policy_type = Policy;
0058
0059    explicit extreme_value_distribution(RealType a = 0, RealType b = 1)
0060       : m_a(a), m_b(b)
0061    {
0062       RealType err;
0063       detail::verify_scale_b("boost::math::extreme_value_distribution<%1%>::extreme_value_distribution", b, &err, Policy());
0064       detail::check_finite("boost::math::extreme_value_distribution<%1%>::extreme_value_distribution", a, &err, Policy());
0065    } // extreme_value_distribution
0066
0067    RealType location()const { return m_a; }
0068    RealType scale()const { return m_b; }
0069
0070 private:
0071    RealType m_a;
0072    RealType m_b;
0073 };
0074
0075 using extreme_value = extreme_value_distribution<double>;
0076
0077 #ifdef __cpp_deduction_guides
0078 template <class RealType>
0079 extreme_value_distribution(RealType)->extreme_value_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0080 template <class RealType>
0081 extreme_value_distribution(RealType,RealType)->extreme_value_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0082 #endif
0083
0084 template <class RealType, class Policy>
0085 inline std::pair<RealType, RealType> range(const extreme_value_distribution<RealType, Policy>& /*dist*/)
0086 { // Range of permissible values for random variable x.
0087    using boost::math::tools::max_value;
0088    return std::pair<RealType, RealType>(
0089       std::numeric_limits<RealType>::has_infinity ? -std::numeric_limits<RealType>::infinity() : -max_value<RealType>(),
0090       std::numeric_limits<RealType>::has_infinity ? std::numeric_limits<RealType>::infinity() : max_value<RealType>());
0091 }
0092
0093 template <class RealType, class Policy>
0094 inline std::pair<RealType, RealType> support(const extreme_value_distribution<RealType, Policy>& /*dist*/)
0095 { // Range of supported values for random variable x.
0096    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0097    using boost::math::tools::max_value;
0098    return std::pair<RealType, RealType>(-max_value<RealType>(),  max_value<RealType>());
0099 }
0100
0101 template <class RealType, class Policy>
0102 inline RealType pdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x)
0103 {
0104    BOOST_MATH_STD_USING // for ADL of std functions
0105
0106    static const char* function = "boost::math::pdf(const extreme_value_distribution<%1%>&, %1%)";
0107
0108    RealType a = dist.location();
0109    RealType b = dist.scale();
0110    RealType result = 0;
0111    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
0112       return result;
0113    if(0 == detail::check_finite(function, a, &result, Policy()))
0114       return result;
0115    if((boost::math::isinf)(x))
0116       return 0.0f;
0117    if(0 == detail::check_x(function, x, &result, Policy()))
0118       return result;
0119    RealType e = (a - x) / b;
0120    if(e < tools::log_max_value<RealType>())
0121       result = exp(e) * exp(-exp(e)) / b;
0122    // else.... result *must* be zero since exp(e) is infinite...
0123    return result;
0124 } // pdf
0125
0126 template <class RealType, class Policy>
0127 inline RealType logpdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x)
0128 {
0129    BOOST_MATH_STD_USING // for ADL of std functions
0130
0131    static const char* function = "boost::math::logpdf(const extreme_value_distribution<%1%>&, %1%)";
0132
0133    RealType a = dist.location();
0134    RealType b = dist.scale();
0135    RealType result = -std::numeric_limits<RealType>::infinity();
0136    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
0137       return result;
0138    if(0 == detail::check_finite(function, a, &result, Policy()))
0139       return result;
0140    if((boost::math::isinf)(x))
0141       return 0.0f;
0142    if(0 == detail::check_x(function, x, &result, Policy()))
0143       return result;
0144    RealType e = (a - x) / b;
0145    if(e < tools::log_max_value<RealType>())
0146       result = log(1/b) + e - exp(e);
0147    // else.... result *must* be zero since exp(e) is infinite...
0148    return result;
0149 } // logpdf
0150
0151 template <class RealType, class Policy>
0152 inline RealType cdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x)
0153 {
0154    BOOST_MATH_STD_USING // for ADL of std functions
0155
0156    static const char* function = "boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)";
0157
0158    if((boost::math::isinf)(x))
0159       return x < 0 ? 0.0f : 1.0f;
0160    RealType a = dist.location();
0161    RealType b = dist.scale();
0162    RealType result = 0;
0163    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
0164       return result;
0165    if(0 == detail::check_finite(function, a, &result, Policy()))
0166       return result;
0167    if(0 == detail::check_finite(function, a, &result, Policy()))
0168       return result;
0169    if(0 == detail::check_x("boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)", x, &result, Policy()))
0170       return result;
0171
0172    result = exp(-exp((a-x)/b));
0173
0174    return result;
0175 } // cdf
0176
0177 template <class RealType, class Policy>
0178 inline RealType logcdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x)
0179 {
0180    BOOST_MATH_STD_USING // for ADL of std functions
0181
0182    static const char* function = "boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)";
0183
0184    if((boost::math::isinf)(x))
0185       return x < 0 ? 0.0f : 1.0f;
0186    RealType a = dist.location();
0187    RealType b = dist.scale();
0188    RealType result = 0;
0189    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
0190       return result;
0191    if(0 == detail::check_finite(function, a, &result, Policy()))
0192       return result;
0193    if(0 == detail::check_finite(function, a, &result, Policy()))
0194       return result;
0195    if(0 == detail::check_x("boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)", x, &result, Policy()))
0196       return result;
0197
0198    result = -exp((a-x)/b);
0199
0200    return result;
0201 } // logcdf
0202
0203 template <class RealType, class Policy>
0204 RealType quantile(const extreme_value_distribution<RealType, Policy>& dist, const RealType& p)
0205 {
0206    BOOST_MATH_STD_USING // for ADL of std functions
0207
0208    static const char* function = "boost::math::quantile(const extreme_value_distribution<%1%>&, %1%)";
0209
0210    RealType a = dist.location();
0211    RealType b = dist.scale();
0212    RealType result = 0;
0213    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
0214       return result;
0215    if(0 == detail::check_finite(function, a, &result, Policy()))
0216       return result;
0217    if(0 == detail::check_probability(function, p, &result, Policy()))
0218       return result;
0219
0220    if(p == 0)
0221       return -policies::raise_overflow_error<RealType>(function, 0, Policy());
0222    if(p == 1)
0223       return policies::raise_overflow_error<RealType>(function, 0, Policy());
0224
0225    result = a - log(-log(p)) * b;
0226
0227    return result;
0228 } // quantile
0229
0230 template <class RealType, class Policy>
0231 inline RealType cdf(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c)
0232 {
0233    BOOST_MATH_STD_USING // for ADL of std functions
0234
0235    static const char* function = "boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)";
0236
0237    if((boost::math::isinf)(c.param))
0238       return c.param < 0 ? 1.0f : 0.0f;
0239    RealType a = c.dist.location();
0240    RealType b = c.dist.scale();
0241    RealType result = 0;
0242    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
0243       return result;
0244    if(0 == detail::check_finite(function, a, &result, Policy()))
0245       return result;
0246    if(0 == detail::check_x(function, c.param, &result, Policy()))
0247       return result;
0248
0249    result = -boost::math::expm1(-exp((a-c.param)/b), Policy());
0250
0251    return result;
0252 }
0253
0254 template <class RealType, class Policy>
0255 inline RealType logcdf(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c)
0256 {
0257    BOOST_MATH_STD_USING // for ADL of std functions
0258
0259    static const char* function = "boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)";
0260
0261    if((boost::math::isinf)(c.param))
0262       return c.param < 0 ? 1.0f : 0.0f;
0263    RealType a = c.dist.location();
0264    RealType b = c.dist.scale();
0265    RealType result = 0;
0266    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
0267       return result;
0268    if(0 == detail::check_finite(function, a, &result, Policy()))
0269       return result;
0270    if(0 == detail::check_x(function, c.param, &result, Policy()))
0271       return result;
0272
0273    result = log1p(-exp(-exp((a-c.param)/b)), Policy());
0274
0275    return result;
0276 }
0277
0278 template <class RealType, class Policy>
0279 RealType quantile(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c)
0280 {
0281    BOOST_MATH_STD_USING // for ADL of std functions
0282
0283    static const char* function = "boost::math::quantile(const extreme_value_distribution<%1%>&, %1%)";
0284
0285    RealType a = c.dist.location();
0286    RealType b = c.dist.scale();
0287    RealType q = c.param;
0288    RealType result = 0;
0289    if(0 == detail::verify_scale_b(function, b, &result, Policy()))
0290       return result;
0291    if(0 == detail::check_finite(function, a, &result, Policy()))
0292       return result;
0293    if(0 == detail::check_probability(function, q, &result, Policy()))
0294       return result;
0295
0296    if(q == 0)
0297       return policies::raise_overflow_error<RealType>(function, 0, Policy());
0298    if(q == 1)
0299       return -policies::raise_overflow_error<RealType>(function, 0, Policy());
0300
0301    result = a - log(-boost::math::log1p(-q, Policy())) * b;
0302
0303    return result;
0304 }
0305
0306 template <class RealType, class Policy>
0307 inline RealType mean(const extreme_value_distribution<RealType, Policy>& dist)
0308 {
0309    RealType a = dist.location();
0310    RealType b = dist.scale();
0311    RealType result = 0;
0312    if(0 == detail::verify_scale_b("boost::math::mean(const extreme_value_distribution<%1%>&)", b, &result, Policy()))
0313       return result;
0314    if (0 == detail::check_finite("boost::math::mean(const extreme_value_distribution<%1%>&)", a, &result, Policy()))
0315       return result;
0316    return a + constants::euler<RealType>() * b;
0317 }
0318
0319 template <class RealType, class Policy>
0320 inline RealType standard_deviation(const extreme_value_distribution<RealType, Policy>& dist)
0321 {
0322    BOOST_MATH_STD_USING // for ADL of std functions.
0323
0324    RealType b = dist.scale();
0325    RealType result = 0;
0326    if(0 == detail::verify_scale_b("boost::math::standard_deviation(const extreme_value_distribution<%1%>&)", b, &result, Policy()))
0327       return result;
0328    if(0 == detail::check_finite("boost::math::standard_deviation(const extreme_value_distribution<%1%>&)", dist.location(), &result, Policy()))
0329       return result;
0330    return constants::pi<RealType>() * b / sqrt(static_cast<RealType>(6));
0331 }
0332
0333 template <class RealType, class Policy>
0334 inline RealType mode(const extreme_value_distribution<RealType, Policy>& dist)
0335 {
0336    return dist.location();
0337 }
0338
0339 template <class RealType, class Policy>
0340 inline RealType median(const extreme_value_distribution<RealType, Policy>& dist)
0341 {
0342   using constants::ln_ln_two;
0343    return dist.location() - dist.scale() * ln_ln_two<RealType>();
0344 }
0345
0346 template <class RealType, class Policy>
0347 inline RealType skewness(const extreme_value_distribution<RealType, Policy>& /*dist*/)
0348 {
0349    //
0350    // This is 12 * sqrt(6) * zeta(3) / pi^3:
0351    // See http://mathworld.wolfram.com/ExtremeValueDistribution.html
0352    //
0353    return static_cast<RealType>(1.1395470994046486574927930193898461120875997958366L);
0354 }
0355
0356 template <class RealType, class Policy>
0357 inline RealType kurtosis(const extreme_value_distribution<RealType, Policy>& /*dist*/)
0358 {
0359    // See http://mathworld.wolfram.com/ExtremeValueDistribution.html
0360    return RealType(27) / 5;
0361 }
0362
0363 template <class RealType, class Policy>
0364 inline RealType kurtosis_excess(const extreme_value_distribution<RealType, Policy>& /*dist*/)
0365 {
0366    // See http://mathworld.wolfram.com/ExtremeValueDistribution.html
0367    return RealType(12) / 5;
0368 }
0369
0370
0371 } // namespace math
0372 } // namespace boost
0373
0374 #ifdef _MSC_VER
0375 # pragma warning(pop)
0376 #endif
0377
0378 // This include must be at the end, *after* the accessors
0379 // for this distribution have been defined, in order to
0380 // keep compilers that support two-phase lookup happy.
0381 #include <boost/math/distributions/detail/derived_accessors.hpp>
0382
0383 #endif // BOOST_STATS_EXTREME_VALUE_HPP