]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | #include <boost/config/no_tr1/cmath.hpp> | |
16 | ||
17 | // | |
18 | // This is the maximum extreme value distribution, see | |
19 | // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366g.htm | |
20 | // and http://mathworld.wolfram.com/ExtremeValueDistribution.html | |
21 | // Also known as a Fisher-Tippett distribution, a log-Weibull | |
22 | // distribution or a Gumbel distribution. | |
23 | ||
24 | #include <utility> | |
25 | ||
26 | #ifdef BOOST_MSVC | |
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 | template <class RealType, class Policy> | |
77 | inline const std::pair<RealType, RealType> range(const extreme_value_distribution<RealType, Policy>& /*dist*/) | |
78 | { // Range of permissible values for random variable x. | |
79 | using boost::math::tools::max_value; | |
80 | return std::pair<RealType, RealType>( | |
81 | std::numeric_limits<RealType>::has_infinity ? -std::numeric_limits<RealType>::infinity() : -max_value<RealType>(), | |
82 | std::numeric_limits<RealType>::has_infinity ? std::numeric_limits<RealType>::infinity() : max_value<RealType>()); | |
83 | } | |
84 | ||
85 | template <class RealType, class Policy> | |
86 | inline const std::pair<RealType, RealType> support(const extreme_value_distribution<RealType, Policy>& /*dist*/) | |
87 | { // Range of supported values for random variable x. | |
88 | // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. | |
89 | using boost::math::tools::max_value; | |
90 | return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); | |
91 | } | |
92 | ||
93 | template <class RealType, class Policy> | |
94 | 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 | ||
98 | static const char* function = "boost::math::pdf(const extreme_value_distribution<%1%>&, %1%)"; | |
99 | ||
100 | RealType a = dist.location(); | |
101 | RealType b = dist.scale(); | |
102 | RealType result = 0; | |
103 | if(0 == detail::verify_scale_b(function, b, &result, Policy())) | |
104 | return result; | |
105 | if(0 == detail::check_finite(function, a, &result, Policy())) | |
106 | return result; | |
107 | if((boost::math::isinf)(x)) | |
108 | return 0.0f; | |
109 | if(0 == detail::check_x(function, x, &result, Policy())) | |
110 | return result; | |
111 | RealType e = (a - x) / b; | |
112 | if(e < tools::log_max_value<RealType>()) | |
113 | result = exp(e) * exp(-exp(e)) / b; | |
114 | // else.... result *must* be zero since exp(e) is infinite... | |
115 | return result; | |
116 | ||
117 | ||
118 | template <class RealType, class Policy> | |
119 | inline RealType cdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x) | |
120 | { | |
121 | BOOST_MATH_STD_USING // for ADL of std functions | |
122 | ||
123 | static const char* function = "boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)"; | |
124 | ||
125 | if((boost::math::isinf)(x)) | |
126 | return x < 0 ? 0.0f : 1.0f; | |
127 | RealType a = dist.location(); | |
128 | RealType b = dist.scale(); | |
129 | RealType result = 0; | |
130 | if(0 == detail::verify_scale_b(function, b, &result, Policy())) | |
131 | return result; | |
132 | if(0 == detail::check_finite(function, a, &result, Policy())) | |
133 | return result; | |
134 | if(0 == detail::check_finite(function, a, &result, Policy())) | |
135 | return result; | |
136 | if(0 == detail::check_x("boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)", x, &result, Policy())) | |
137 | return result; | |
138 | ||
139 | result = exp(-exp((a-x)/b)); | |
140 | ||
141 | return result; | |
142 | } // cdf | |
143 | ||
144 | template <class RealType, class Policy> | |
145 | RealType quantile(const extreme_value_distribution<RealType, Policy>& dist, const RealType& p) | |
146 | { | |
147 | BOOST_MATH_STD_USING // for ADL of std functions | |
148 | ||
149 | static const char* function = "boost::math::quantile(const extreme_value_distribution<%1%>&, %1%)"; | |
150 | ||
151 | RealType a = dist.location(); | |
152 | RealType b = dist.scale(); | |
153 | RealType result = 0; | |
154 | if(0 == detail::verify_scale_b(function, b, &result, Policy())) | |
155 | return result; | |
156 | if(0 == detail::check_finite(function, a, &result, Policy())) | |
157 | return result; | |
158 | if(0 == detail::check_probability(function, p, &result, Policy())) | |
159 | return result; | |
160 | ||
161 | if(p == 0) | |
162 | return -policies::raise_overflow_error<RealType>(function, 0, Policy()); | |
163 | if(p == 1) | |
164 | return policies::raise_overflow_error<RealType>(function, 0, Policy()); | |
165 | ||
166 | result = a - log(-log(p)) * b; | |
167 | ||
168 | return result; | |
169 | } // quantile | |
170 | ||
171 | template <class RealType, class Policy> | |
172 | inline RealType cdf(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c) | |
173 | { | |
174 | BOOST_MATH_STD_USING // for ADL of std functions | |
175 | ||
176 | static const char* function = "boost::math::cdf(const extreme_value_distribution<%1%>&, %1%)"; | |
177 | ||
178 | if((boost::math::isinf)(c.param)) | |
179 | return c.param < 0 ? 1.0f : 0.0f; | |
180 | RealType a = c.dist.location(); | |
181 | RealType b = c.dist.scale(); | |
182 | RealType result = 0; | |
183 | if(0 == detail::verify_scale_b(function, b, &result, Policy())) | |
184 | return result; | |
185 | if(0 == detail::check_finite(function, a, &result, Policy())) | |
186 | return result; | |
187 | if(0 == detail::check_x(function, c.param, &result, Policy())) | |
188 | return result; | |
189 | ||
190 | result = -boost::math::expm1(-exp((a-c.param)/b), Policy()); | |
191 | ||
192 | return result; | |
193 | } | |
194 | ||
195 | template <class RealType, class Policy> | |
196 | RealType quantile(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c) | |
197 | { | |
198 | BOOST_MATH_STD_USING // for ADL of std functions | |
199 | ||
200 | static const char* function = "boost::math::quantile(const extreme_value_distribution<%1%>&, %1%)"; | |
201 | ||
202 | RealType a = c.dist.location(); | |
203 | RealType b = c.dist.scale(); | |
204 | RealType q = c.param; | |
205 | RealType result = 0; | |
206 | if(0 == detail::verify_scale_b(function, b, &result, Policy())) | |
207 | return result; | |
208 | if(0 == detail::check_finite(function, a, &result, Policy())) | |
209 | return result; | |
210 | if(0 == detail::check_probability(function, q, &result, Policy())) | |
211 | return result; | |
212 | ||
213 | if(q == 0) | |
214 | return policies::raise_overflow_error<RealType>(function, 0, Policy()); | |
215 | if(q == 1) | |
216 | return -policies::raise_overflow_error<RealType>(function, 0, Policy()); | |
217 | ||
218 | result = a - log(-boost::math::log1p(-q, Policy())) * b; | |
219 | ||
220 | return result; | |
221 | } | |
222 | ||
223 | template <class RealType, class Policy> | |
224 | inline RealType mean(const extreme_value_distribution<RealType, Policy>& dist) | |
225 | { | |
226 | RealType a = dist.location(); | |
227 | RealType b = dist.scale(); | |
228 | RealType result = 0; | |
229 | if(0 == detail::verify_scale_b("boost::math::mean(const extreme_value_distribution<%1%>&)", b, &result, Policy())) | |
230 | return result; | |
231 | if(0 == detail::check_scale("boost::math::mean(const extreme_value_distribution<%1%>&)", a, &result, Policy())) | |
232 | return result; | |
233 | return a + constants::euler<RealType>() * b; | |
234 | } | |
235 | ||
236 | template <class RealType, class Policy> | |
237 | inline RealType standard_deviation(const extreme_value_distribution<RealType, Policy>& dist) | |
238 | { | |
239 | BOOST_MATH_STD_USING // for ADL of std functions. | |
240 | ||
241 | RealType b = dist.scale(); | |
242 | RealType result = 0; | |
243 | if(0 == detail::verify_scale_b("boost::math::standard_deviation(const extreme_value_distribution<%1%>&)", b, &result, Policy())) | |
244 | return result; | |
245 | if(0 == detail::check_scale("boost::math::standard_deviation(const extreme_value_distribution<%1%>&)", dist.location(), &result, Policy())) | |
246 | return result; | |
247 | return constants::pi<RealType>() * b / sqrt(static_cast<RealType>(6)); | |
248 | } | |
249 | ||
250 | template <class RealType, class Policy> | |
251 | inline RealType mode(const extreme_value_distribution<RealType, Policy>& dist) | |
252 | { | |
253 | return dist.location(); | |
254 | } | |
255 | ||
256 | template <class RealType, class Policy> | |
257 | inline RealType median(const extreme_value_distribution<RealType, Policy>& dist) | |
258 | { | |
259 | using constants::ln_ln_two; | |
260 | return dist.location() - dist.scale() * ln_ln_two<RealType>(); | |
261 | } | |
262 | ||
263 | template <class RealType, class Policy> | |
264 | inline RealType skewness(const extreme_value_distribution<RealType, Policy>& /*dist*/) | |
265 | { | |
266 | // | |
267 | // This is 12 * sqrt(6) * zeta(3) / pi^3: | |
268 | // See http://mathworld.wolfram.com/ExtremeValueDistribution.html | |
269 | // | |
270 | return static_cast<RealType>(1.1395470994046486574927930193898461120875997958366L); | |
271 | } | |
272 | ||
273 | template <class RealType, class Policy> | |
274 | inline RealType kurtosis(const extreme_value_distribution<RealType, Policy>& /*dist*/) | |
275 | { | |
276 | // See http://mathworld.wolfram.com/ExtremeValueDistribution.html | |
277 | return RealType(27) / 5; | |
278 | } | |
279 | ||
280 | template <class RealType, class Policy> | |
281 | inline RealType kurtosis_excess(const extreme_value_distribution<RealType, Policy>& /*dist*/) | |
282 | { | |
283 | // See http://mathworld.wolfram.com/ExtremeValueDistribution.html | |
284 | return RealType(12) / 5; | |
285 | } | |
286 | ||
287 | ||
288 | } // namespace math | |
289 | } // namespace boost | |
290 | ||
291 | #ifdef BOOST_MSVC | |
292 | # pragma warning(pop) | |
293 | #endif | |
294 | ||
295 | // This include must be at the end, *after* the accessors | |
296 | // for this distribution have been defined, in order to | |
297 | // keep compilers that support two-phase lookup happy. | |
298 | #include <boost/math/distributions/detail/derived_accessors.hpp> | |
299 | ||
300 | #endif // BOOST_STATS_EXTREME_VALUE_HPP |