]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright John Maddock 2007. |
2 | // Copyright Paul A. Bristow 2007, 2009 | |
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 | ||
7 | #ifndef BOOST_STATS_PARETO_HPP | |
8 | #define BOOST_STATS_PARETO_HPP | |
9 | ||
10 | // http://en.wikipedia.org/wiki/Pareto_distribution | |
11 | // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm | |
12 | // Also: | |
13 | // Weisstein, Eric W. "Pareto Distribution." | |
14 | // From MathWorld--A Wolfram Web Resource. | |
15 | // http://mathworld.wolfram.com/ParetoDistribution.html | |
16 | // Handbook of Statistical Distributions with Applications, K Krishnamoorthy, ISBN 1-58488-635-8, Chapter 23, pp 257 - 267. | |
17 | // Caution KK's a and b are the reverse of Mathworld! | |
18 | ||
19 | #include <boost/math/distributions/fwd.hpp> | |
20 | #include <boost/math/distributions/complement.hpp> | |
21 | #include <boost/math/distributions/detail/common_error_handling.hpp> | |
22 | #include <boost/math/special_functions/powm1.hpp> | |
23 | ||
24 | #include <utility> // for BOOST_CURRENT_VALUE? | |
25 | ||
26 | namespace boost | |
27 | { | |
28 | namespace math | |
29 | { | |
30 | namespace detail | |
31 | { // Parameter checking. | |
32 | template <class RealType, class Policy> | |
33 | inline bool check_pareto_scale( | |
34 | const char* function, | |
35 | RealType scale, | |
36 | RealType* result, const Policy& pol) | |
37 | { | |
38 | if((boost::math::isfinite)(scale)) | |
39 | { // any > 0 finite value is OK. | |
40 | if (scale > 0) | |
41 | { | |
42 | return true; | |
43 | } | |
44 | else | |
45 | { | |
46 | *result = policies::raise_domain_error<RealType>( | |
47 | function, | |
48 | "Scale parameter is %1%, but must be > 0!", scale, pol); | |
49 | return false; | |
50 | } | |
51 | } | |
52 | else | |
53 | { // Not finite. | |
54 | *result = policies::raise_domain_error<RealType>( | |
55 | function, | |
56 | "Scale parameter is %1%, but must be finite!", scale, pol); | |
57 | return false; | |
58 | } | |
59 | } // bool check_pareto_scale | |
60 | ||
61 | template <class RealType, class Policy> | |
62 | inline bool check_pareto_shape( | |
63 | const char* function, | |
64 | RealType shape, | |
65 | RealType* result, const Policy& pol) | |
66 | { | |
67 | if((boost::math::isfinite)(shape)) | |
68 | { // Any finite value > 0 is OK. | |
69 | if (shape > 0) | |
70 | { | |
71 | return true; | |
72 | } | |
73 | else | |
74 | { | |
75 | *result = policies::raise_domain_error<RealType>( | |
76 | function, | |
77 | "Shape parameter is %1%, but must be > 0!", shape, pol); | |
78 | return false; | |
79 | } | |
80 | } | |
81 | else | |
82 | { // Not finite. | |
83 | *result = policies::raise_domain_error<RealType>( | |
84 | function, | |
85 | "Shape parameter is %1%, but must be finite!", shape, pol); | |
86 | return false; | |
87 | } | |
88 | } // bool check_pareto_shape( | |
89 | ||
90 | template <class RealType, class Policy> | |
91 | inline bool check_pareto_x( | |
92 | const char* function, | |
93 | RealType const& x, | |
94 | RealType* result, const Policy& pol) | |
95 | { | |
96 | if((boost::math::isfinite)(x)) | |
97 | { // | |
98 | if (x > 0) | |
99 | { | |
100 | return true; | |
101 | } | |
102 | else | |
103 | { | |
104 | *result = policies::raise_domain_error<RealType>( | |
105 | function, | |
106 | "x parameter is %1%, but must be > 0 !", x, pol); | |
107 | return false; | |
108 | } | |
109 | } | |
110 | else | |
111 | { // Not finite.. | |
112 | *result = policies::raise_domain_error<RealType>( | |
113 | function, | |
114 | "x parameter is %1%, but must be finite!", x, pol); | |
115 | return false; | |
116 | } | |
117 | } // bool check_pareto_x | |
118 | ||
119 | template <class RealType, class Policy> | |
120 | inline bool check_pareto( // distribution parameters. | |
121 | const char* function, | |
122 | RealType scale, | |
123 | RealType shape, | |
124 | RealType* result, const Policy& pol) | |
125 | { | |
126 | return check_pareto_scale(function, scale, result, pol) | |
127 | && check_pareto_shape(function, shape, result, pol); | |
128 | } // bool check_pareto( | |
129 | ||
130 | } // namespace detail | |
131 | ||
132 | template <class RealType = double, class Policy = policies::policy<> > | |
133 | class pareto_distribution | |
134 | { | |
135 | public: | |
136 | typedef RealType value_type; | |
137 | typedef Policy policy_type; | |
138 | ||
139 | pareto_distribution(RealType l_scale = 1, RealType l_shape = 1) | |
140 | : m_scale(l_scale), m_shape(l_shape) | |
141 | { // Constructor. | |
142 | RealType result = 0; | |
143 | detail::check_pareto("boost::math::pareto_distribution<%1%>::pareto_distribution", l_scale, l_shape, &result, Policy()); | |
144 | } | |
145 | ||
146 | RealType scale()const | |
147 | { // AKA Xm and Wolfram b and beta | |
148 | return m_scale; | |
149 | } | |
150 | ||
151 | RealType shape()const | |
152 | { // AKA k and Wolfram a and alpha | |
153 | return m_shape; | |
154 | } | |
155 | private: | |
156 | // Data members: | |
157 | RealType m_scale; // distribution scale (xm) or beta | |
158 | RealType m_shape; // distribution shape (k) or alpha | |
159 | }; | |
160 | ||
161 | typedef pareto_distribution<double> pareto; // Convenience to allow pareto(2., 3.); | |
162 | ||
163 | template <class RealType, class Policy> | |
164 | inline const std::pair<RealType, RealType> range(const pareto_distribution<RealType, Policy>& /*dist*/) | |
165 | { // Range of permissible values for random variable x. | |
166 | using boost::math::tools::max_value; | |
167 | return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // scale zero to + infinity. | |
168 | } // range | |
169 | ||
170 | template <class RealType, class Policy> | |
171 | inline const std::pair<RealType, RealType> support(const pareto_distribution<RealType, Policy>& dist) | |
172 | { // Range of supported values for random variable x. | |
173 | // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. | |
174 | using boost::math::tools::max_value; | |
175 | return std::pair<RealType, RealType>(dist.scale(), max_value<RealType>() ); // scale to + infinity. | |
176 | } // support | |
177 | ||
178 | template <class RealType, class Policy> | |
179 | inline RealType pdf(const pareto_distribution<RealType, Policy>& dist, const RealType& x) | |
180 | { | |
181 | BOOST_MATH_STD_USING // for ADL of std function pow. | |
182 | static const char* function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)"; | |
183 | RealType scale = dist.scale(); | |
184 | RealType shape = dist.shape(); | |
185 | RealType result = 0; | |
186 | if(false == (detail::check_pareto_x(function, x, &result, Policy()) | |
187 | && detail::check_pareto(function, scale, shape, &result, Policy()))) | |
188 | return result; | |
189 | if (x < scale) | |
190 | { // regardless of shape, pdf is zero (or should be disallow x < scale and throw an exception?). | |
191 | return 0; | |
192 | } | |
193 | result = shape * pow(scale, shape) / pow(x, shape+1); | |
194 | return result; | |
195 | ||
196 | ||
197 | template <class RealType, class Policy> | |
198 | inline RealType cdf(const pareto_distribution<RealType, Policy>& dist, const RealType& x) | |
199 | { | |
200 | BOOST_MATH_STD_USING // for ADL of std function pow. | |
201 | static const char* function = "boost::math::cdf(const pareto_distribution<%1%>&, %1%)"; | |
202 | RealType scale = dist.scale(); | |
203 | RealType shape = dist.shape(); | |
204 | RealType result = 0; | |
205 | ||
206 | if(false == (detail::check_pareto_x(function, x, &result, Policy()) | |
207 | && detail::check_pareto(function, scale, shape, &result, Policy()))) | |
208 | return result; | |
209 | ||
210 | if (x <= scale) | |
211 | { // regardless of shape, cdf is zero. | |
212 | return 0; | |
213 | } | |
214 | ||
215 | // result = RealType(1) - pow((scale / x), shape); | |
216 | result = -boost::math::powm1(scale/x, shape, Policy()); // should be more accurate. | |
217 | return result; | |
218 | } // cdf | |
219 | ||
220 | template <class RealType, class Policy> | |
221 | inline RealType quantile(const pareto_distribution<RealType, Policy>& dist, const RealType& p) | |
222 | { | |
223 | BOOST_MATH_STD_USING // for ADL of std function pow. | |
224 | static const char* function = "boost::math::quantile(const pareto_distribution<%1%>&, %1%)"; | |
225 | RealType result = 0; | |
226 | RealType scale = dist.scale(); | |
227 | RealType shape = dist.shape(); | |
228 | if(false == (detail::check_probability(function, p, &result, Policy()) | |
229 | && detail::check_pareto(function, scale, shape, &result, Policy()))) | |
230 | { | |
231 | return result; | |
232 | } | |
233 | if (p == 0) | |
234 | { | |
235 | return scale; // x must be scale (or less). | |
236 | } | |
237 | if (p == 1) | |
238 | { | |
239 | return policies::raise_overflow_error<RealType>(function, 0, Policy()); // x = + infinity. | |
240 | } | |
241 | result = scale / | |
242 | (pow((1 - p), 1 / shape)); | |
243 | // K. Krishnamoorthy, ISBN 1-58488-635-8 eq 23.1.3 | |
244 | return result; | |
245 | } // quantile | |
246 | ||
247 | template <class RealType, class Policy> | |
248 | inline RealType cdf(const complemented2_type<pareto_distribution<RealType, Policy>, RealType>& c) | |
249 | { | |
250 | BOOST_MATH_STD_USING // for ADL of std function pow. | |
251 | static const char* function = "boost::math::cdf(const pareto_distribution<%1%>&, %1%)"; | |
252 | RealType result = 0; | |
253 | RealType x = c.param; | |
254 | RealType scale = c.dist.scale(); | |
255 | RealType shape = c.dist.shape(); | |
256 | if(false == (detail::check_pareto_x(function, x, &result, Policy()) | |
257 | && detail::check_pareto(function, scale, shape, &result, Policy()))) | |
258 | return result; | |
259 | ||
260 | if (x <= scale) | |
261 | { // regardless of shape, cdf is zero, and complement is unity. | |
262 | return 1; | |
263 | } | |
264 | result = pow((scale/x), shape); | |
265 | ||
266 | return result; | |
267 | } // cdf complement | |
268 | ||
269 | template <class RealType, class Policy> | |
270 | inline RealType quantile(const complemented2_type<pareto_distribution<RealType, Policy>, RealType>& c) | |
271 | { | |
272 | BOOST_MATH_STD_USING // for ADL of std function pow. | |
273 | static const char* function = "boost::math::quantile(const pareto_distribution<%1%>&, %1%)"; | |
274 | RealType result = 0; | |
275 | RealType q = c.param; | |
276 | RealType scale = c.dist.scale(); | |
277 | RealType shape = c.dist.shape(); | |
278 | if(false == (detail::check_probability(function, q, &result, Policy()) | |
279 | && detail::check_pareto(function, scale, shape, &result, Policy()))) | |
280 | { | |
281 | return result; | |
282 | } | |
283 | if (q == 1) | |
284 | { | |
285 | return scale; // x must be scale (or less). | |
286 | } | |
287 | if (q == 0) | |
288 | { | |
289 | return policies::raise_overflow_error<RealType>(function, 0, Policy()); // x = + infinity. | |
290 | } | |
291 | result = scale / (pow(q, 1 / shape)); | |
292 | // K. Krishnamoorthy, ISBN 1-58488-635-8 eq 23.1.3 | |
293 | return result; | |
294 | } // quantile complement | |
295 | ||
296 | template <class RealType, class Policy> | |
297 | inline RealType mean(const pareto_distribution<RealType, Policy>& dist) | |
298 | { | |
299 | RealType result = 0; | |
300 | static const char* function = "boost::math::mean(const pareto_distribution<%1%>&, %1%)"; | |
301 | if(false == detail::check_pareto(function, dist.scale(), dist.shape(), &result, Policy())) | |
302 | { | |
303 | return result; | |
304 | } | |
305 | if (dist.shape() > RealType(1)) | |
306 | { | |
307 | return dist.shape() * dist.scale() / (dist.shape() - 1); | |
308 | } | |
309 | else | |
310 | { | |
311 | using boost::math::tools::max_value; | |
312 | return max_value<RealType>(); // +infinity. | |
313 | } | |
314 | } // mean | |
315 | ||
316 | template <class RealType, class Policy> | |
317 | inline RealType mode(const pareto_distribution<RealType, Policy>& dist) | |
318 | { | |
319 | return dist.scale(); | |
320 | } // mode | |
321 | ||
322 | template <class RealType, class Policy> | |
323 | inline RealType median(const pareto_distribution<RealType, Policy>& dist) | |
324 | { | |
325 | RealType result = 0; | |
326 | static const char* function = "boost::math::median(const pareto_distribution<%1%>&, %1%)"; | |
327 | if(false == detail::check_pareto(function, dist.scale(), dist.shape(), &result, Policy())) | |
328 | { | |
329 | return result; | |
330 | } | |
331 | BOOST_MATH_STD_USING | |
332 | return dist.scale() * pow(RealType(2), (1/dist.shape())); | |
333 | } // median | |
334 | ||
335 | template <class RealType, class Policy> | |
336 | inline RealType variance(const pareto_distribution<RealType, Policy>& dist) | |
337 | { | |
338 | RealType result = 0; | |
339 | RealType scale = dist.scale(); | |
340 | RealType shape = dist.shape(); | |
341 | static const char* function = "boost::math::variance(const pareto_distribution<%1%>&, %1%)"; | |
342 | if(false == detail::check_pareto(function, scale, shape, &result, Policy())) | |
343 | { | |
344 | return result; | |
345 | } | |
346 | if (shape > 2) | |
347 | { | |
348 | result = (scale * scale * shape) / | |
349 | ((shape - 1) * (shape - 1) * (shape - 2)); | |
350 | } | |
351 | else | |
352 | { | |
353 | result = policies::raise_domain_error<RealType>( | |
354 | function, | |
355 | "variance is undefined for shape <= 2, but got %1%.", dist.shape(), Policy()); | |
356 | } | |
357 | return result; | |
358 | } // variance | |
359 | ||
360 | template <class RealType, class Policy> | |
361 | inline RealType skewness(const pareto_distribution<RealType, Policy>& dist) | |
362 | { | |
363 | BOOST_MATH_STD_USING | |
364 | RealType result = 0; | |
365 | RealType shape = dist.shape(); | |
366 | static const char* function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)"; | |
367 | if(false == detail::check_pareto(function, dist.scale(), shape, &result, Policy())) | |
368 | { | |
369 | return result; | |
370 | } | |
371 | if (shape > 3) | |
372 | { | |
373 | result = sqrt((shape - 2) / shape) * | |
374 | 2 * (shape + 1) / | |
375 | (shape - 3); | |
376 | } | |
377 | else | |
378 | { | |
379 | result = policies::raise_domain_error<RealType>( | |
380 | function, | |
381 | "skewness is undefined for shape <= 3, but got %1%.", dist.shape(), Policy()); | |
382 | } | |
383 | return result; | |
384 | } // skewness | |
385 | ||
386 | template <class RealType, class Policy> | |
387 | inline RealType kurtosis(const pareto_distribution<RealType, Policy>& dist) | |
388 | { | |
389 | RealType result = 0; | |
390 | RealType shape = dist.shape(); | |
391 | static const char* function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)"; | |
392 | if(false == detail::check_pareto(function, dist.scale(), shape, &result, Policy())) | |
393 | { | |
394 | return result; | |
395 | } | |
396 | if (shape > 4) | |
397 | { | |
398 | result = 3 * ((shape - 2) * (3 * shape * shape + shape + 2)) / | |
399 | (shape * (shape - 3) * (shape - 4)); | |
400 | } | |
401 | else | |
402 | { | |
403 | result = policies::raise_domain_error<RealType>( | |
404 | function, | |
405 | "kurtosis_excess is undefined for shape <= 4, but got %1%.", shape, Policy()); | |
406 | } | |
407 | return result; | |
408 | } // kurtosis | |
409 | ||
410 | template <class RealType, class Policy> | |
411 | inline RealType kurtosis_excess(const pareto_distribution<RealType, Policy>& dist) | |
412 | { | |
413 | RealType result = 0; | |
414 | RealType shape = dist.shape(); | |
415 | static const char* function = "boost::math::pdf(const pareto_distribution<%1%>&, %1%)"; | |
416 | if(false == detail::check_pareto(function, dist.scale(), shape, &result, Policy())) | |
417 | { | |
418 | return result; | |
419 | } | |
420 | if (shape > 4) | |
421 | { | |
422 | result = 6 * ((shape * shape * shape) + (shape * shape) - 6 * shape - 2) / | |
423 | (shape * (shape - 3) * (shape - 4)); | |
424 | } | |
425 | else | |
426 | { | |
427 | result = policies::raise_domain_error<RealType>( | |
428 | function, | |
429 | "kurtosis_excess is undefined for shape <= 4, but got %1%.", dist.shape(), Policy()); | |
430 | } | |
431 | return result; | |
432 | } // kurtosis_excess | |
433 | ||
434 | } // namespace math | |
435 | } // namespace boost | |
436 | ||
437 | // This include must be at the end, *after* the accessors | |
438 | // for this distribution have been defined, in order to | |
439 | // keep compilers that support two-phase lookup happy. | |
440 | #include <boost/math/distributions/detail/derived_accessors.hpp> | |
441 | ||
442 | #endif // BOOST_STATS_PARETO_HPP | |
443 | ||
444 |