393 lines · plain
1// Copyright John Maddock 2006.2// Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_STATS_EXTREME_VALUE_HPP8#define BOOST_STATS_EXTREME_VALUE_HPP9 10#include <boost/math/tools/config.hpp>11#include <boost/math/tools/numeric_limits.hpp>12#include <boost/math/tools/tuple.hpp>13#include <boost/math/tools/precision.hpp>14#include <boost/math/constants/constants.hpp>15#include <boost/math/special_functions/log1p.hpp>16#include <boost/math/special_functions/expm1.hpp>17#include <boost/math/distributions/complement.hpp>18#include <boost/math/distributions/detail/common_error_handling.hpp>19#include <boost/math/policies/policy.hpp>20#include <boost/math/policies/error_handling.hpp>21 22//23// This is the maximum extreme value distribution, see24// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366g.htm25// and http://mathworld.wolfram.com/ExtremeValueDistribution.html26// Also known as a Fisher-Tippett distribution, a log-Weibull27// distribution or a Gumbel distribution.28 29#ifndef BOOST_MATH_HAS_NVRTC30#include <boost/math/distributions/fwd.hpp>31#include <utility>32#include <cmath>33#endif34 35#ifdef _MSC_VER36# pragma warning(push)37# pragma warning(disable: 4702) // unreachable code (return after domain_error throw).38#endif39 40namespace boost{ namespace math{41 42namespace detail{43//44// Error check:45//46template <class RealType, class Policy>47BOOST_MATH_GPU_ENABLED inline bool verify_scale_b(const char* function, RealType b, RealType* presult, const Policy& pol)48{49 if((b <= 0) || !(boost::math::isfinite)(b))50 {51 *presult = policies::raise_domain_error<RealType>(52 function,53 "The scale parameter \"b\" must be finite and > 0, but was: %1%.", b, pol);54 return false;55 }56 return true;57}58 59} // namespace detail60 61template <class RealType = double, class Policy = policies::policy<> >62class extreme_value_distribution63{64public:65 using value_type = RealType;66 using policy_type = Policy;67 68 BOOST_MATH_GPU_ENABLED explicit extreme_value_distribution(RealType a = 0, RealType b = 1)69 : m_a(a), m_b(b)70 {71 RealType err;72 detail::verify_scale_b("boost::math::extreme_value_distribution<%1%>::extreme_value_distribution", b, &err, Policy());73 detail::check_finite("boost::math::extreme_value_distribution<%1%>::extreme_value_distribution", a, &err, Policy());74 } // extreme_value_distribution75 76 BOOST_MATH_GPU_ENABLED RealType location()const { return m_a; }77 BOOST_MATH_GPU_ENABLED RealType scale()const { return m_b; }78 79private:80 RealType m_a;81 RealType m_b;82};83 84using extreme_value = extreme_value_distribution<double>;85 86#ifdef __cpp_deduction_guides87template <class RealType>88extreme_value_distribution(RealType)->extreme_value_distribution<typename boost::math::tools::promote_args<RealType>::type>;89template <class RealType>90extreme_value_distribution(RealType,RealType)->extreme_value_distribution<typename boost::math::tools::promote_args<RealType>::type>;91#endif92 93template <class RealType, class Policy>94BOOST_MATH_GPU_ENABLED inline boost::math::pair<RealType, RealType> range(const extreme_value_distribution<RealType, Policy>& /*dist*/)95{ // Range of permissible values for random variable x.96 using boost::math::tools::max_value;97 return boost::math::pair<RealType, RealType>(98 boost::math::numeric_limits<RealType>::has_infinity ? -boost::math::numeric_limits<RealType>::infinity() : -max_value<RealType>(), 99 boost::math::numeric_limits<RealType>::has_infinity ? boost::math::numeric_limits<RealType>::infinity() : max_value<RealType>());100}101 102template <class RealType, class Policy>103BOOST_MATH_GPU_ENABLED inline boost::math::pair<RealType, RealType> support(const extreme_value_distribution<RealType, Policy>& /*dist*/)104{ // Range of supported values for random variable x.105 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.106 using boost::math::tools::max_value;107 return boost::math::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());108}109 110template <class RealType, class Policy>111BOOST_MATH_GPU_ENABLED inline RealType pdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x)112{113 BOOST_MATH_STD_USING // for ADL of std functions114 115 constexpr auto function = "boost::math::pdf(const extreme_value_distribution<%1%>&, %1%)";116 117 RealType a = dist.location();118 RealType b = dist.scale();119 RealType result = 0;120 if(0 == detail::verify_scale_b(function, b, &result, Policy()))121 return result;122 if(0 == detail::check_finite(function, a, &result, Policy()))123 return result;124 if((boost::math::isinf)(x))125 return 0.0f;126 if(0 == detail::check_x(function, x, &result, Policy()))127 return result;128 RealType e = (a - x) / b;129 if(e < tools::log_max_value<RealType>())130 result = exp(e) * exp(-exp(e)) / b;131 // else.... result *must* be zero since exp(e) is infinite...132 return result;133} // pdf134 135template <class RealType, class Policy>136BOOST_MATH_GPU_ENABLED inline RealType logpdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x)137{138 BOOST_MATH_STD_USING // for ADL of std functions139 140 constexpr auto function = "boost::math::logpdf(const extreme_value_distribution<%1%>&, %1%)";141 142 RealType a = dist.location();143 RealType b = dist.scale();144 RealType result = -boost::math::numeric_limits<RealType>::infinity();145 if(0 == detail::verify_scale_b(function, b, &result, Policy()))146 return result;147 if(0 == detail::check_finite(function, a, &result, Policy()))148 return result;149 if((boost::math::isinf)(x))150 return 0.0f;151 if(0 == detail::check_x(function, x, &result, Policy()))152 return result;153 RealType e = (a - x) / b;154 if(e < tools::log_max_value<RealType>())155 result = log(1/b) + e - exp(e);156 // else.... result *must* be zero since exp(e) is infinite...157 return result;158} // logpdf159 160template <class RealType, class Policy>161BOOST_MATH_GPU_ENABLED inline RealType cdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x)162{163 BOOST_MATH_STD_USING // for ADL of std functions164 165 constexpr auto function = "boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)";166 167 if((boost::math::isinf)(x))168 return x < 0 ? 0.0f : 1.0f;169 RealType a = dist.location();170 RealType b = dist.scale();171 RealType result = 0;172 if(0 == detail::verify_scale_b(function, b, &result, Policy()))173 return result;174 if(0 == detail::check_finite(function, a, &result, Policy()))175 return result;176 if(0 == detail::check_finite(function, a, &result, Policy()))177 return result;178 if(0 == detail::check_x("boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)", x, &result, Policy()))179 return result;180 181 result = exp(-exp((a-x)/b));182 183 return result;184} // cdf185 186template <class RealType, class Policy>187BOOST_MATH_GPU_ENABLED inline RealType logcdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x)188{189 BOOST_MATH_STD_USING // for ADL of std functions190 191 constexpr auto function = "boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)";192 193 if((boost::math::isinf)(x))194 return x < 0 ? 0.0f : 1.0f;195 RealType a = dist.location();196 RealType b = dist.scale();197 RealType result = 0;198 if(0 == detail::verify_scale_b(function, b, &result, Policy()))199 return result;200 if(0 == detail::check_finite(function, a, &result, Policy()))201 return result;202 if(0 == detail::check_finite(function, a, &result, Policy()))203 return result;204 if(0 == detail::check_x("boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)", x, &result, Policy()))205 return result;206 207 result = -exp((a-x)/b);208 209 return result;210} // logcdf211 212template <class RealType, class Policy>213BOOST_MATH_GPU_ENABLED RealType quantile(const extreme_value_distribution<RealType, Policy>& dist, const RealType& p)214{215 BOOST_MATH_STD_USING // for ADL of std functions216 217 constexpr auto function = "boost::math::quantile(const extreme_value_distribution<%1%>&, %1%)";218 219 RealType a = dist.location();220 RealType b = dist.scale();221 RealType result = 0;222 if(0 == detail::verify_scale_b(function, b, &result, Policy()))223 return result;224 if(0 == detail::check_finite(function, a, &result, Policy()))225 return result;226 if(0 == detail::check_probability(function, p, &result, Policy()))227 return result;228 229 if(p == 0)230 return -policies::raise_overflow_error<RealType>(function, 0, Policy());231 if(p == 1)232 return policies::raise_overflow_error<RealType>(function, 0, Policy());233 234 result = a - log(-log(p)) * b;235 236 return result;237} // quantile238 239template <class RealType, class Policy>240BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c)241{242 BOOST_MATH_STD_USING // for ADL of std functions243 244 constexpr auto function = "boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)";245 246 if((boost::math::isinf)(c.param))247 return c.param < 0 ? 1.0f : 0.0f;248 RealType a = c.dist.location();249 RealType b = c.dist.scale();250 RealType result = 0;251 if(0 == detail::verify_scale_b(function, b, &result, Policy()))252 return result;253 if(0 == detail::check_finite(function, a, &result, Policy()))254 return result;255 if(0 == detail::check_x(function, c.param, &result, Policy()))256 return result;257 258 result = -boost::math::expm1(-exp((a-c.param)/b), Policy());259 260 return result;261}262 263template <class RealType, class Policy>264BOOST_MATH_GPU_ENABLED inline RealType logcdf(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c)265{266 BOOST_MATH_STD_USING // for ADL of std functions267 268 constexpr auto function = "boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)";269 270 if((boost::math::isinf)(c.param))271 return c.param < 0 ? 1.0f : 0.0f;272 RealType a = c.dist.location();273 RealType b = c.dist.scale();274 RealType result = 0;275 if(0 == detail::verify_scale_b(function, b, &result, Policy()))276 return result;277 if(0 == detail::check_finite(function, a, &result, Policy()))278 return result;279 if(0 == detail::check_x(function, c.param, &result, Policy()))280 return result;281 282 result = log1p(-exp(-exp((a-c.param)/b)), Policy());283 284 return result;285}286 287template <class RealType, class Policy>288BOOST_MATH_GPU_ENABLED RealType quantile(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c)289{290 BOOST_MATH_STD_USING // for ADL of std functions291 292 constexpr auto function = "boost::math::quantile(const extreme_value_distribution<%1%>&, %1%)";293 294 RealType a = c.dist.location();295 RealType b = c.dist.scale();296 RealType q = c.param;297 RealType result = 0;298 if(0 == detail::verify_scale_b(function, b, &result, Policy()))299 return result;300 if(0 == detail::check_finite(function, a, &result, Policy()))301 return result;302 if(0 == detail::check_probability(function, q, &result, Policy()))303 return result;304 305 if(q == 0)306 return policies::raise_overflow_error<RealType>(function, 0, Policy());307 if(q == 1)308 return -policies::raise_overflow_error<RealType>(function, 0, Policy());309 310 result = a - log(-boost::math::log1p(-q, Policy())) * b;311 312 return result;313}314 315template <class RealType, class Policy>316BOOST_MATH_GPU_ENABLED inline RealType mean(const extreme_value_distribution<RealType, Policy>& dist)317{318 RealType a = dist.location();319 RealType b = dist.scale();320 RealType result = 0;321 if(0 == detail::verify_scale_b("boost::math::mean(const extreme_value_distribution<%1%>&)", b, &result, Policy()))322 return result;323 if (0 == detail::check_finite("boost::math::mean(const extreme_value_distribution<%1%>&)", a, &result, Policy()))324 return result;325 return a + constants::euler<RealType>() * b;326}327 328template <class RealType, class Policy>329BOOST_MATH_GPU_ENABLED inline RealType standard_deviation(const extreme_value_distribution<RealType, Policy>& dist)330{331 BOOST_MATH_STD_USING // for ADL of std functions.332 333 RealType b = dist.scale();334 RealType result = 0;335 if(0 == detail::verify_scale_b("boost::math::standard_deviation(const extreme_value_distribution<%1%>&)", b, &result, Policy()))336 return result;337 if(0 == detail::check_finite("boost::math::standard_deviation(const extreme_value_distribution<%1%>&)", dist.location(), &result, Policy()))338 return result;339 return constants::pi<RealType>() * b / sqrt(static_cast<RealType>(6));340}341 342template <class RealType, class Policy>343BOOST_MATH_GPU_ENABLED inline RealType mode(const extreme_value_distribution<RealType, Policy>& dist)344{345 return dist.location();346}347 348template <class RealType, class Policy>349BOOST_MATH_GPU_ENABLED inline RealType median(const extreme_value_distribution<RealType, Policy>& dist)350{351 using constants::ln_ln_two;352 return dist.location() - dist.scale() * ln_ln_two<RealType>();353}354 355template <class RealType, class Policy>356BOOST_MATH_GPU_ENABLED inline RealType skewness(const extreme_value_distribution<RealType, Policy>& /*dist*/)357{358 //359 // This is 12 * sqrt(6) * zeta(3) / pi^3:360 // See http://mathworld.wolfram.com/ExtremeValueDistribution.html361 //362 return static_cast<RealType>(1.1395470994046486574927930193898461120875997958366L);363}364 365template <class RealType, class Policy>366BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const extreme_value_distribution<RealType, Policy>& /*dist*/)367{368 // See http://mathworld.wolfram.com/ExtremeValueDistribution.html369 return RealType(27) / 5;370}371 372template <class RealType, class Policy>373BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const extreme_value_distribution<RealType, Policy>& /*dist*/)374{375 // See http://mathworld.wolfram.com/ExtremeValueDistribution.html376 return RealType(12) / 5;377}378 379 380} // namespace math381} // namespace boost382 383#ifdef _MSC_VER384# pragma warning(pop)385#endif386 387// This include must be at the end, *after* the accessors388// for this distribution have been defined, in order to389// keep compilers that support two-phase lookup happy.390#include <boost/math/distributions/detail/derived_accessors.hpp>391 392#endif // BOOST_STATS_EXTREME_VALUE_HPP393