extreme_value.hpp 14 KB

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