| 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 | using value_type = RealType; |
| 57 | using policy_type = Policy; |
| 58 | |
| 59 | explicit 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; |
| 72 | RealType m_b; |
| 73 | }; |
| 74 | |
| 75 | using extreme_value = extreme_value_distribution<double>; |
| 76 | |
| 77 | #ifdef __cpp_deduction_guides |
| 78 | template <class RealType> |
| 79 | extreme_value_distribution(RealType)->extreme_value_distribution<typename boost::math::tools::promote_args<RealType>::type>; |
| 80 | template <class RealType> |
| 81 | extreme_value_distribution(RealType,RealType)->extreme_value_distribution<typename boost::math::tools::promote_args<RealType>::type>; |
| 82 | #endif |
| 83 | |
| 84 | template <class RealType, class Policy> |
| 85 | inline std::pair<RealType, RealType> range(const extreme_value_distribution<RealType, Policy>& /*dist*/) |
| 86 | { // Range of permissible values for random variable x. |
| 87 | using boost::math::tools::max_value; |
| 88 | return std::pair<RealType, RealType>( |
| 89 | std::numeric_limits<RealType>::has_infinity ? -std::numeric_limits<RealType>::infinity() : -max_value<RealType>(), |
| 90 | std::numeric_limits<RealType>::has_infinity ? std::numeric_limits<RealType>::infinity() : max_value<RealType>()); |
| 91 | } |
| 92 | |
| 93 | template <class RealType, class Policy> |
| 94 | inline std::pair<RealType, RealType> support(const extreme_value_distribution<RealType, Policy>& /*dist*/) |
| 95 | { // Range of supported values for random variable x. |
| 96 | // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. |
| 97 | using boost::math::tools::max_value; |
| 98 | return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); |
| 99 | } |
| 100 | |
| 101 | template <class RealType, class Policy> |
| 102 | inline RealType pdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x) |
| 103 | { |
| 104 | BOOST_MATH_STD_USING // for ADL of std functions |
| 105 | |
| 106 | static const char* function = "boost::math::pdf(const extreme_value_distribution<%1%>&, %1%)" ; |
| 107 | |
| 108 | RealType a = dist.location(); |
| 109 | RealType b = dist.scale(); |
| 110 | RealType result = 0; |
| 111 | if(0 == detail::verify_scale_b(function, b, &result, Policy())) |
| 112 | return result; |
| 113 | if(0 == detail::check_finite(function, a, &result, Policy())) |
| 114 | return result; |
| 115 | if((boost::math::isinf)(x)) |
| 116 | return 0.0f; |
| 117 | if(0 == detail::check_x(function, x, &result, Policy())) |
| 118 | return result; |
| 119 | RealType e = (a - x) / b; |
| 120 | if(e < tools::log_max_value<RealType>()) |
| 121 | result = exp(e) * exp(-exp(e)) / b; |
| 122 | // else.... result *must* be zero since exp(e) is infinite... |
| 123 | return result; |
| 124 | } // pdf |
| 125 | |
| 126 | template <class RealType, class Policy> |
| 127 | inline RealType logpdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x) |
| 128 | { |
| 129 | BOOST_MATH_STD_USING // for ADL of std functions |
| 130 | |
| 131 | static const char* function = "boost::math::logpdf(const extreme_value_distribution<%1%>&, %1%)" ; |
| 132 | |
| 133 | RealType a = dist.location(); |
| 134 | RealType b = dist.scale(); |
| 135 | RealType result = -std::numeric_limits<RealType>::infinity(); |
| 136 | if(0 == detail::verify_scale_b(function, b, &result, Policy())) |
| 137 | return result; |
| 138 | if(0 == detail::check_finite(function, a, &result, Policy())) |
| 139 | return result; |
| 140 | if((boost::math::isinf)(x)) |
| 141 | return 0.0f; |
| 142 | if(0 == detail::check_x(function, x, &result, Policy())) |
| 143 | return result; |
| 144 | RealType e = (a - x) / b; |
| 145 | if(e < tools::log_max_value<RealType>()) |
| 146 | result = log(1/b) + e - exp(e); |
| 147 | // else.... result *must* be zero since exp(e) is infinite... |
| 148 | return result; |
| 149 | } // logpdf |
| 150 | |
| 151 | template <class RealType, class Policy> |
| 152 | inline RealType cdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x) |
| 153 | { |
| 154 | BOOST_MATH_STD_USING // for ADL of std functions |
| 155 | |
| 156 | static const char* function = "boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)" ; |
| 157 | |
| 158 | if((boost::math::isinf)(x)) |
| 159 | return x < 0 ? 0.0f : 1.0f; |
| 160 | RealType a = dist.location(); |
| 161 | RealType b = dist.scale(); |
| 162 | RealType result = 0; |
| 163 | if(0 == detail::verify_scale_b(function, b, &result, Policy())) |
| 164 | return result; |
| 165 | if(0 == detail::check_finite(function, a, &result, Policy())) |
| 166 | return result; |
| 167 | if(0 == detail::check_finite(function, a, &result, Policy())) |
| 168 | return result; |
| 169 | if(0 == detail::check_x("boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)" , x, &result, Policy())) |
| 170 | return result; |
| 171 | |
| 172 | result = exp(-exp((a-x)/b)); |
| 173 | |
| 174 | return result; |
| 175 | } // cdf |
| 176 | |
| 177 | template <class RealType, class Policy> |
| 178 | inline RealType logcdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x) |
| 179 | { |
| 180 | BOOST_MATH_STD_USING // for ADL of std functions |
| 181 | |
| 182 | static const char* function = "boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)" ; |
| 183 | |
| 184 | if((boost::math::isinf)(x)) |
| 185 | return x < 0 ? 0.0f : 1.0f; |
| 186 | RealType a = dist.location(); |
| 187 | RealType b = dist.scale(); |
| 188 | RealType result = 0; |
| 189 | if(0 == detail::verify_scale_b(function, b, &result, Policy())) |
| 190 | return result; |
| 191 | if(0 == detail::check_finite(function, a, &result, Policy())) |
| 192 | return result; |
| 193 | if(0 == detail::check_finite(function, a, &result, Policy())) |
| 194 | return result; |
| 195 | if(0 == detail::check_x("boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)" , x, &result, Policy())) |
| 196 | return result; |
| 197 | |
| 198 | result = -exp((a-x)/b); |
| 199 | |
| 200 | return result; |
| 201 | } // logcdf |
| 202 | |
| 203 | template <class RealType, class Policy> |
| 204 | RealType quantile(const extreme_value_distribution<RealType, Policy>& dist, const RealType& p) |
| 205 | { |
| 206 | BOOST_MATH_STD_USING // for ADL of std functions |
| 207 | |
| 208 | static const char* function = "boost::math::quantile(const extreme_value_distribution<%1%>&, %1%)" ; |
| 209 | |
| 210 | RealType a = dist.location(); |
| 211 | RealType b = dist.scale(); |
| 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, p, &result, Policy())) |
| 218 | return result; |
| 219 | |
| 220 | if(p == 0) |
| 221 | return -policies::raise_overflow_error<RealType>(function, 0, Policy()); |
| 222 | if(p == 1) |
| 223 | return policies::raise_overflow_error<RealType>(function, 0, Policy()); |
| 224 | |
| 225 | result = a - log(-log(p)) * b; |
| 226 | |
| 227 | return result; |
| 228 | } // quantile |
| 229 | |
| 230 | template <class RealType, class Policy> |
| 231 | inline RealType cdf(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c) |
| 232 | { |
| 233 | BOOST_MATH_STD_USING // for ADL of std functions |
| 234 | |
| 235 | static const char* function = "boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)" ; |
| 236 | |
| 237 | if((boost::math::isinf)(c.param)) |
| 238 | return c.param < 0 ? 1.0f : 0.0f; |
| 239 | RealType a = c.dist.location(); |
| 240 | RealType b = c.dist.scale(); |
| 241 | RealType result = 0; |
| 242 | if(0 == detail::verify_scale_b(function, b, &result, Policy())) |
| 243 | return result; |
| 244 | if(0 == detail::check_finite(function, a, &result, Policy())) |
| 245 | return result; |
| 246 | if(0 == detail::check_x(function, c.param, &result, Policy())) |
| 247 | return result; |
| 248 | |
| 249 | result = -boost::math::expm1(-exp((a-c.param)/b), Policy()); |
| 250 | |
| 251 | return result; |
| 252 | } |
| 253 | |
| 254 | template <class RealType, class Policy> |
| 255 | inline RealType logcdf(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c) |
| 256 | { |
| 257 | BOOST_MATH_STD_USING // for ADL of std functions |
| 258 | |
| 259 | static const char* function = "boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)" ; |
| 260 | |
| 261 | if((boost::math::isinf)(c.param)) |
| 262 | return c.param < 0 ? 1.0f : 0.0f; |
| 263 | RealType a = c.dist.location(); |
| 264 | RealType b = c.dist.scale(); |
| 265 | RealType result = 0; |
| 266 | if(0 == detail::verify_scale_b(function, b, &result, Policy())) |
| 267 | return result; |
| 268 | if(0 == detail::check_finite(function, a, &result, Policy())) |
| 269 | return result; |
| 270 | if(0 == detail::check_x(function, c.param, &result, Policy())) |
| 271 | return result; |
| 272 | |
| 273 | result = log1p(-exp(-exp((a-c.param)/b)), Policy()); |
| 274 | |
| 275 | return result; |
| 276 | } |
| 277 | |
| 278 | template <class RealType, class Policy> |
| 279 | RealType quantile(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c) |
| 280 | { |
| 281 | BOOST_MATH_STD_USING // for ADL of std functions |
| 282 | |
| 283 | static const char* function = "boost::math::quantile(const extreme_value_distribution<%1%>&, %1%)" ; |
| 284 | |
| 285 | RealType a = c.dist.location(); |
| 286 | RealType b = c.dist.scale(); |
| 287 | RealType q = c.param; |
| 288 | RealType result = 0; |
| 289 | if(0 == detail::verify_scale_b(function, b, &result, Policy())) |
| 290 | return result; |
| 291 | if(0 == detail::check_finite(function, a, &result, Policy())) |
| 292 | return result; |
| 293 | if(0 == detail::check_probability(function, q, &result, Policy())) |
| 294 | return result; |
| 295 | |
| 296 | if(q == 0) |
| 297 | return policies::raise_overflow_error<RealType>(function, 0, Policy()); |
| 298 | if(q == 1) |
| 299 | return -policies::raise_overflow_error<RealType>(function, 0, Policy()); |
| 300 | |
| 301 | result = a - log(-boost::math::log1p(-q, Policy())) * b; |
| 302 | |
| 303 | return result; |
| 304 | } |
| 305 | |
| 306 | template <class RealType, class Policy> |
| 307 | inline RealType mean(const extreme_value_distribution<RealType, Policy>& dist) |
| 308 | { |
| 309 | RealType a = dist.location(); |
| 310 | RealType b = dist.scale(); |
| 311 | RealType result = 0; |
| 312 | if(0 == detail::verify_scale_b("boost::math::mean(const extreme_value_distribution<%1%>&)" , b, &result, Policy())) |
| 313 | return result; |
| 314 | if (0 == detail::check_finite("boost::math::mean(const extreme_value_distribution<%1%>&)" , a, &result, Policy())) |
| 315 | return result; |
| 316 | return a + constants::euler<RealType>() * b; |
| 317 | } |
| 318 | |
| 319 | template <class RealType, class Policy> |
| 320 | inline RealType standard_deviation(const extreme_value_distribution<RealType, Policy>& dist) |
| 321 | { |
| 322 | BOOST_MATH_STD_USING // for ADL of std functions. |
| 323 | |
| 324 | RealType b = dist.scale(); |
| 325 | RealType result = 0; |
| 326 | if(0 == detail::verify_scale_b("boost::math::standard_deviation(const extreme_value_distribution<%1%>&)" , b, &result, Policy())) |
| 327 | return result; |
| 328 | if(0 == detail::check_finite("boost::math::standard_deviation(const extreme_value_distribution<%1%>&)" , dist.location(), &result, Policy())) |
| 329 | return result; |
| 330 | return constants::pi<RealType>() * b / sqrt(static_cast<RealType>(6)); |
| 331 | } |
| 332 | |
| 333 | template <class RealType, class Policy> |
| 334 | inline RealType mode(const extreme_value_distribution<RealType, Policy>& dist) |
| 335 | { |
| 336 | return dist.location(); |
| 337 | } |
| 338 | |
| 339 | template <class RealType, class Policy> |
| 340 | inline RealType median(const extreme_value_distribution<RealType, Policy>& dist) |
| 341 | { |
| 342 | using constants::ln_ln_two; |
| 343 | return dist.location() - dist.scale() * ln_ln_two<RealType>(); |
| 344 | } |
| 345 | |
| 346 | template <class RealType, class Policy> |
| 347 | inline RealType skewness(const extreme_value_distribution<RealType, Policy>& /*dist*/) |
| 348 | { |
| 349 | // |
| 350 | // This is 12 * sqrt(6) * zeta(3) / pi^3: |
| 351 | // See http://mathworld.wolfram.com/ExtremeValueDistribution.html |
| 352 | // |
| 353 | return static_cast<RealType>(1.1395470994046486574927930193898461120875997958366L); |
| 354 | } |
| 355 | |
| 356 | template <class RealType, class Policy> |
| 357 | inline RealType kurtosis(const extreme_value_distribution<RealType, Policy>& /*dist*/) |
| 358 | { |
| 359 | // See http://mathworld.wolfram.com/ExtremeValueDistribution.html |
| 360 | return RealType(27) / 5; |
| 361 | } |
| 362 | |
| 363 | template <class RealType, class Policy> |
| 364 | inline RealType kurtosis_excess(const extreme_value_distribution<RealType, Policy>& /*dist*/) |
| 365 | { |
| 366 | // See http://mathworld.wolfram.com/ExtremeValueDistribution.html |
| 367 | return RealType(12) / 5; |
| 368 | } |
| 369 | |
| 370 | |
| 371 | } // namespace math |
| 372 | } // namespace boost |
| 373 | |
| 374 | #ifdef _MSC_VER |
| 375 | # pragma warning(pop) |
| 376 | #endif |
| 377 | |
| 378 | // This include must be at the end, *after* the accessors |
| 379 | // for this distribution have been defined, in order to |
| 380 | // keep compilers that support two-phase lookup happy. |
| 381 | #include <boost/math/distributions/detail/derived_accessors.hpp> |
| 382 | |
| 383 | #endif // BOOST_STATS_EXTREME_VALUE_HPP |
| 384 | |