1 // Copyright Paul A. Bristow 2007, 2009.
2 // Copyright John Maddock 2006.
3 // Use, modification and distribution are subject to the
4 // Boost Software License, Version 1.0.
5 // (See accompanying file LICENSE_1_0.txt
6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
10 // http://en.wikipedia.org/wiki/pareto_distribution
11 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm
13 // Weisstein, Eric W. "pareto Distribution."
14 // From MathWorld--A Wolfram Web Resource.
15 // http://mathworld.wolfram.com/paretoDistribution.html
19 # pragma warning(disable: 4127) // conditional expression is constant.
20 # pragma warning (disable : 4996) // POSIX name for this item is deprecated
21 # pragma warning (disable : 4224) // nonstandard extension used : formal parameter 'arg' was previously defined as a type
22 # pragma warning (disable : 4180) // qualifier applied to function type has no meaning; ignored
23 # pragma warning(disable: 4100) // unreferenced formal parameter.
26 #include <boost/math/tools/test.hpp> // for real_concept
27 #include <boost/math/concepts/real_concept.hpp> // for real_concept
28 #define BOOST_TEST_MAIN
29 #include <boost/test/unit_test.hpp> // Boost.Test
30 #include <boost/test/tools/floating_point_comparison.hpp>
32 #include <boost/math/distributions/pareto.hpp>
33 using boost::math::pareto_distribution
;
34 #include <boost/math/tools/test.hpp>
35 #include "test_out_of_range.hpp"
40 using std::setprecision
;
42 using std::numeric_limits
;
44 template <class RealType
>
45 void check_pareto(RealType scale
, RealType shape
, RealType x
, RealType p
, RealType q
, RealType tol
)
47 BOOST_CHECK_CLOSE_FRACTION(
49 pareto_distribution
<RealType
>(scale
, shape
), // distribution.
50 x
), // random variable.
52 tol
); // tolerance eps.
53 BOOST_CHECK_CLOSE_FRACTION(
56 pareto_distribution
<RealType
>(scale
, shape
), // distribution.
57 x
)), // random variable.
58 q
, // probability complement.
59 tol
); // tolerance eps.
60 BOOST_CHECK_CLOSE_FRACTION(
61 ::boost::math::quantile(
62 pareto_distribution
<RealType
>(scale
, shape
), // distribution.
64 x
, // random variable.
65 tol
); // tolerance eps.
66 BOOST_CHECK_CLOSE_FRACTION(
67 ::boost::math::quantile(
69 pareto_distribution
<RealType
>(scale
, shape
), // distribution.
70 q
)), // probability complement.
71 x
, // random variable.
72 tol
); // tolerance eps.
75 template <class RealType
>
76 void test_spots(RealType
)
78 // Basic sanity checks.
80 // Tolerance are based on units of epsilon, but capped at
81 // double precision, since that's the limit of our test data:
83 RealType tol
= (std::max
)((RealType
)boost::math::tools::epsilon
<double>(), boost::math::tools::epsilon
<RealType
>());
84 RealType tol5eps
= tol
* 5;
85 RealType tol10eps
= tol
* 10;
86 RealType tol100eps
= tol
* 100;
87 RealType tol1000eps
= tol
* 1000;
90 static_cast<RealType
>(1.1L), //
91 static_cast<RealType
>(5.5L),
92 static_cast<RealType
>(2.2L),
93 static_cast<RealType
>(0.97790291308792L),
94 static_cast<RealType
>(0.0220970869120796L),
98 static_cast<RealType
>(0.5L),
99 static_cast<RealType
>(10.1L),
100 static_cast<RealType
>(1.5L),
101 static_cast<RealType
>(0.99998482686481L),
102 static_cast<RealType
>(1.51731351900608e-005L),
103 tol100eps
* 1000); // Much less accurate as p close to unity.
106 static_cast<RealType
>(0.1L),
107 static_cast<RealType
>(2.3L),
108 static_cast<RealType
>(1.5L),
109 static_cast<RealType
>(0.99802762220697L),
110 static_cast<RealType
>(0.00197237779302972L),
113 // Example from 23.3 page 259
115 static_cast<RealType
>(2.30444301457005L),
116 static_cast<RealType
>(4),
117 static_cast<RealType
>(2.4L),
118 static_cast<RealType
>(0.15L),
119 static_cast<RealType
>(0.85L),
123 static_cast<RealType
>(2),
124 static_cast<RealType
>(3),
125 static_cast<RealType
>(3.4L),
126 static_cast<RealType
>(0.796458375737838L),
127 static_cast<RealType
>(0.203541624262162L),
130 check_pareto( // Probability near 0.5
131 static_cast<RealType
>(2),
132 static_cast<RealType
>(2),
133 static_cast<RealType
>(3),
134 static_cast<RealType
>(0.5555555555555555555555555555555555555556L),
135 static_cast<RealType
>(0.4444444444444444444444444444444444444444L),
136 tol5eps
); // accurate.
141 // pdf for shapes 1, 2 & 3 (exact)
142 BOOST_CHECK_CLOSE_FRACTION(
143 pdf(pareto_distribution
<RealType
>(1, 1), 1),
144 static_cast<RealType
>(1), //
147 BOOST_CHECK_CLOSE_FRACTION( pdf(pareto_distribution
<RealType
>(1, 2), 1),
148 static_cast<RealType
>(2), //
151 BOOST_CHECK_CLOSE_FRACTION( pdf(pareto_distribution
<RealType
>(1, 3), 1),
152 static_cast<RealType
>(3), //
156 BOOST_CHECK_EQUAL( // x = scale
157 cdf(pareto_distribution
<RealType
>(1, 1), 1),
158 static_cast<RealType
>(0) );
160 // Compare with values from StatCalc K. Krishnamoorthy, ISBN 1-58488-635-8 eq 23.1.3
161 BOOST_CHECK_CLOSE_FRACTION( // small x
162 cdf(pareto_distribution
<RealType
>(2, 5), static_cast<RealType
>(3.4)),
163 static_cast<RealType
>(0.929570372227626L), tol5eps
);
165 BOOST_CHECK_CLOSE_FRACTION( // small x
166 cdf(pareto_distribution
<RealType
>(2, 5), static_cast<RealType
>(3.4)),
167 static_cast<RealType
>(1 - 0.0704296277723743L), tol5eps
);
169 BOOST_CHECK_CLOSE_FRACTION( // small x
170 cdf(complement(pareto_distribution
<RealType
>(2, 5), static_cast<RealType
>(3.4))),
171 static_cast<RealType
>(0.0704296277723743L), tol5eps
);
174 BOOST_CHECK_EQUAL( // x = scale
175 quantile(pareto_distribution
<RealType
>(1, 1), 0),
176 static_cast<RealType
>(1) );
178 BOOST_CHECK_EQUAL( // x = scale
179 quantile(complement(pareto_distribution
<RealType
>(1, 1), 1)),
180 static_cast<RealType
>(1) );
182 BOOST_CHECK_CLOSE_FRACTION( // small x
183 cdf(complement(pareto_distribution
<RealType
>(2, 5), static_cast<RealType
>(3.4))),
184 static_cast<RealType
>(0.0704296277723743L), tol5eps
);
186 using namespace std
; // ADL of std names.
188 pareto_distribution
<RealType
> pareto15(1, 5);
189 // Note: shape must be big enough (5) that all moments up to kurtosis are defined
190 // to allow all functions to be tested.
193 BOOST_CHECK_CLOSE_FRACTION(
194 mean(pareto15
), static_cast<RealType
>(1.25), tol5eps
); // 1.25 == 5/4
196 mean(pareto15
), static_cast<RealType
>(1.25)); // 1.25 == 5/4 (expect exact so check equal)
198 pareto_distribution
<RealType
> p12(1, 2); //
200 mean(p12
), static_cast<RealType
>(2)); // Exactly two.
203 BOOST_CHECK_CLOSE_FRACTION(
204 variance(pareto15
), static_cast<RealType
>(0.10416666666666667L), tol5eps
);
206 BOOST_CHECK_CLOSE_FRACTION(
207 standard_deviation(pareto15
), static_cast<RealType
>(0.32274861218395140L), tol5eps
);
208 // hazard: No independent test values found yet.
209 //BOOST_CHECK_CLOSE_FRACTION(
210 // hazard(pareto15, x), pdf(pareto15, x) / cdf(complement(pareto15, x)), tol5eps);
211 //// cumulative hazard:
212 //BOOST_CHECK_CLOSE_FRACTION(
213 // chf(pareto15, x), -log(cdf(complement(pareto15, x))), tol5eps);
214 //// coefficient_of_variation:
215 BOOST_CHECK_CLOSE_FRACTION(
216 coefficient_of_variation(pareto15
), static_cast<RealType
>(0.25819888974716110L), tol5eps
);
218 BOOST_CHECK_CLOSE_FRACTION(
219 mode(pareto15
), static_cast<RealType
>(1), tol5eps
);
221 BOOST_CHECK_CLOSE_FRACTION(
222 median(pareto15
), static_cast<RealType
>(1.1486983549970351L), tol5eps
);
225 BOOST_CHECK_CLOSE_FRACTION(
226 skewness(pareto15
), static_cast<RealType
>(4.6475800154489004L), tol5eps
);
228 BOOST_CHECK_CLOSE_FRACTION(
229 kurtosis(pareto15
), static_cast<RealType
>(73.8L), tol5eps
);
231 BOOST_CHECK_CLOSE_FRACTION(
232 kurtosis_excess(pareto15
), static_cast<RealType
>(70.8L), tol5eps
);
233 // Check difference between kurtosis and excess:
234 BOOST_CHECK_CLOSE_FRACTION(
235 kurtosis_excess(pareto15
), kurtosis(pareto15
) - static_cast<RealType
>(3L), tol5eps
);
236 // Check kurtosis excess = kurtosis - 3;
238 RealType expected_entropy
= 1 + RealType(1)/RealType(5) + log(RealType(1)/RealType(5));
239 BOOST_CHECK_CLOSE_FRACTION(
240 entropy(pareto15
), expected_entropy
, tol5eps
);
242 // Error condition checks:
243 check_out_of_range
<pareto_distribution
<RealType
> >(1, 1);
244 BOOST_MATH_CHECK_THROW(pdf(pareto_distribution
<RealType
>(0, 1), 0), std::domain_error
);
245 BOOST_MATH_CHECK_THROW(pdf(pareto_distribution
<RealType
>(1, 0), 0), std::domain_error
);
246 BOOST_MATH_CHECK_THROW(pdf(pareto_distribution
<RealType
>(-1, 1), 0), std::domain_error
);
247 BOOST_MATH_CHECK_THROW(pdf(pareto_distribution
<RealType
>(1, -1), 0), std::domain_error
);
249 BOOST_MATH_CHECK_THROW(pdf(pareto_distribution
<RealType
>(1, 1), 0), std::domain_error
);
250 BOOST_MATH_CHECK_THROW(cdf(pareto_distribution
<RealType
>(1, 1), 0), std::domain_error
);
252 BOOST_MATH_CHECK_THROW(quantile(pareto_distribution
<RealType
>(1, 1), -1), std::domain_error
);
253 BOOST_MATH_CHECK_THROW(quantile(pareto_distribution
<RealType
>(1, 1), 2), std::domain_error
);
254 } // template <class RealType>void test_spots(RealType)
256 BOOST_AUTO_TEST_CASE( test_main
)
258 // Check that can generate pareto distribution using the two convenience methods:
259 boost::math::pareto
myp1(1., 1); // Using typedef
260 pareto_distribution
<> myp2(1., 1); // Using default RealType double.
261 boost::math::pareto pareto11
; // Use default values (scale = 1, shape = 1).
262 // Note NOT pareto11() as the compiler will interpret as a function!
263 // Basic sanity-check spot values.
265 BOOST_CHECK_EQUAL(pareto11
.scale(), 1); // Check defaults again.
266 BOOST_CHECK_EQUAL(pareto11
.shape(), 1);
267 BOOST_CHECK_EQUAL(myp1
.scale(), 1);
268 BOOST_CHECK_EQUAL(myp1
.shape(), 1);
269 BOOST_CHECK_EQUAL(myp2
.scale(), 1);
270 BOOST_CHECK_EQUAL(myp2
.shape(), 1);
272 // Test range and support using double only,
273 // because it supports numeric_limits max for pseudo-infinity.
274 BOOST_CHECK_EQUAL(range(myp2
).first
, 0); // range 0 to +infinity
275 BOOST_CHECK_EQUAL(range(myp2
).second
, (numeric_limits
<double>::max
)());
276 BOOST_CHECK_EQUAL(support(myp2
).first
, myp2
.scale()); // support scale to + infinity.
277 BOOST_CHECK_EQUAL(support(myp2
).second
, (numeric_limits
<double>::max
)());
279 // Check some bad parameters to the distribution.
280 #ifndef BOOST_NO_EXCEPTIONS
281 BOOST_MATH_CHECK_THROW(boost::math::pareto
mypm1(-1, 1), std::domain_error
); // Using typedef
282 BOOST_MATH_CHECK_THROW(boost::math::pareto
myp0(0, 1), std::domain_error
); // Using typedef
283 BOOST_MATH_CHECK_THROW(boost::math::pareto
myp1m1(1, -1), std::domain_error
); // Using typedef
284 BOOST_MATH_CHECK_THROW(boost::math::pareto
myp10(1, 0), std::domain_error
); // Using typedef
286 BOOST_MATH_CHECK_THROW(boost::math::pareto(-1, 1), std::domain_error
); // Using typedef
287 BOOST_MATH_CHECK_THROW(boost::math::pareto(0, 1), std::domain_error
); // Using typedef
288 BOOST_MATH_CHECK_THROW(boost::math::pareto(1, -1), std::domain_error
); // Using typedef
289 BOOST_MATH_CHECK_THROW(boost::math::pareto(1, 0), std::domain_error
); // Using typedef
292 // Check some moments that should fail because shape not big enough.
293 BOOST_MATH_CHECK_THROW(variance(myp2
), std::domain_error
);
294 BOOST_MATH_CHECK_THROW(standard_deviation(myp2
), std::domain_error
);
295 BOOST_MATH_CHECK_THROW(skewness(myp2
), std::domain_error
);
296 BOOST_MATH_CHECK_THROW(kurtosis(myp2
), std::domain_error
);
297 BOOST_MATH_CHECK_THROW(kurtosis_excess(myp2
), std::domain_error
);
299 // Test on extreme values of distribution parameters,
300 // using just double because it has numeric_limit infinity etc.
301 #ifndef BOOST_NO_EXCEPTIONS
302 BOOST_MATH_CHECK_THROW(boost::math::pareto
mypinf1(+std::numeric_limits
<double>::infinity(), 1), std::domain_error
); // Using typedef
303 BOOST_MATH_CHECK_THROW(boost::math::pareto
myp1inf(1, +std::numeric_limits
<double>::infinity()), std::domain_error
); // Using typedef
304 BOOST_MATH_CHECK_THROW(boost::math::pareto
mypinf1(+std::numeric_limits
<double>::infinity(), +std::numeric_limits
<double>::infinity()), std::domain_error
); // Using typedef
306 BOOST_MATH_CHECK_THROW(boost::math::pareto(+std::numeric_limits
<double>::infinity(), 1), std::domain_error
); // Using typedef
307 BOOST_MATH_CHECK_THROW(boost::math::pareto(1, +std::numeric_limits
<double>::infinity()), std::domain_error
); // Using typedef
308 BOOST_MATH_CHECK_THROW(boost::math::pareto(+std::numeric_limits
<double>::infinity(), +std::numeric_limits
<double>::infinity()), std::domain_error
); // Using typedef
311 // Test on extreme values of random variate x, using just double because it has numeric_limit infinity etc..
312 // No longer allow x to be + or - infinity, then these tests should throw.
313 BOOST_MATH_CHECK_THROW(pdf(pareto11
, +std::numeric_limits
<double>::infinity()), std::domain_error
); // x = + infinity
314 BOOST_MATH_CHECK_THROW(pdf(pareto11
, -std::numeric_limits
<double>::infinity()), std::domain_error
); // x = - infinity
315 BOOST_MATH_CHECK_THROW(cdf(pareto11
, +std::numeric_limits
<double>::infinity()), std::domain_error
); // x = + infinity
316 BOOST_MATH_CHECK_THROW(cdf(pareto11
, -std::numeric_limits
<double>::infinity()), std::domain_error
); // x = - infinity
318 BOOST_CHECK_EQUAL(pdf(pareto11
, 0.5), 0); // x < scale but > 0
319 BOOST_CHECK_EQUAL(pdf(pareto11
, (std::numeric_limits
<double>::min
)()), 0); // x almost zero but > 0
320 BOOST_CHECK_EQUAL(pdf(pareto11
, 1), 1); // x == scale, result == shape == 1
321 BOOST_CHECK_EQUAL(pdf(pareto11
, +(std::numeric_limits
<double>::max
)()), 0); // x = +max, pdf has fallen to zero.
323 BOOST_MATH_CHECK_THROW(pdf(pareto11
, 0), std::domain_error
); // x == 0
324 BOOST_MATH_CHECK_THROW(pdf(pareto11
, -1), std::domain_error
); // x = -1
325 BOOST_MATH_CHECK_THROW(pdf(pareto11
, -(std::numeric_limits
<double>::max
)()), std::domain_error
); // x = - max
326 BOOST_MATH_CHECK_THROW(pdf(pareto11
, -(std::numeric_limits
<double>::min
)()), std::domain_error
); // x = - min
328 BOOST_CHECK_EQUAL(cdf(pareto11
, 1), 0); // x == scale, cdf = zero.
329 BOOST_CHECK_EQUAL(cdf(pareto11
, +(std::numeric_limits
<double>::max
)()), 1); // x = + max, cdf = unity.
331 BOOST_MATH_CHECK_THROW(cdf(pareto11
, 0), std::domain_error
); // x == 0
332 BOOST_MATH_CHECK_THROW(cdf(pareto11
, -(std::numeric_limits
<double>::min
)()), std::domain_error
); // x = - min,
333 BOOST_MATH_CHECK_THROW(cdf(pareto11
, -(std::numeric_limits
<double>::max
)()), std::domain_error
); // x = - max,
335 // (Parameter value, arbitrarily zero, only communicates the floating point type).
336 test_spots(0.0F
); // Test float. OK at decdigits = 0 tol5eps = 0.0001 %
337 test_spots(0.0); // Test double. OK at decdigits 7, tol5eps = 1e07 %
338 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
339 test_spots(0.0L); // Test long double.
340 #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x0582))
341 test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
344 std::cout
<< "<note>The long double tests have been disabled on this platform "
345 "either because the long double overloads of the usual math functions are "
346 "not available at all, or because they are too inaccurate for these tests "
347 "to pass.</note>" << std::endl
;
351 } // BOOST_AUTO_TEST_CASE( test_main )
360 Embedding manifest...
361 Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_pareto.exe"
362 Running 1 test case...
363 *** No errors detected