1 // Copyright Paul A. Bristow 2012.
2 // Copyright John Maddock 2012.
3 // Copyright Benjamin Sobotta 2012
5 // Use, modification and distribution are subject to the
6 // Boost Software License, Version 1.0.
7 // (See accompanying file LICENSE_1_0.txt
8 // or copy at http://www.boost.org/LICENSE_1_0.txt)
11 # pragma warning (disable : 4127) // conditional expression is constant.
12 # pragma warning (disable : 4305) // 'initializing' : truncation from 'double' to 'const float'.
13 # pragma warning (disable : 4310) // cast truncates constant value.
14 # pragma warning (disable : 4512) // assignment operator could not be generated.
17 //#include <pch.hpp> // include directory libs/math/src/tr1/ is needed.
19 #include <boost/math/concepts/real_concept.hpp> // for real_concept
20 #include <boost/test/test_exec_monitor.hpp> // Boost.Test
21 #include <boost/test/floating_point_comparison.hpp>
23 #include <boost/math/distributions/skew_normal.hpp>
24 using boost::math::skew_normal_distribution;
25 using boost::math::skew_normal;
30 using std::setprecision;
32 using std::numeric_limits;
34 template <class RealType>
35 void check_skew_normal(RealType mean, RealType scale, RealType shape, RealType x, RealType p, RealType q, RealType tol)
37 using boost::math::skew_normal_distribution;
39 BOOST_CHECK_CLOSE_FRACTION(
40 ::boost::math::cdf( // Check cdf
41 skew_normal_distribution<RealType>(mean, scale, shape), // distribution.
42 x), // random variable.
45 BOOST_CHECK_CLOSE_FRACTION(
46 ::boost::math::cdf( // Check cdf complement
48 skew_normal_distribution<RealType>(mean, scale, shape), // distribution.
49 x)), // random variable.
50 q, // probability complement.
52 BOOST_CHECK_CLOSE_FRACTION(
53 ::boost::math::quantile( // Check quantile
54 skew_normal_distribution<RealType>(mean, scale, shape), // distribution.
56 x, // random variable.
58 BOOST_CHECK_CLOSE_FRACTION(
59 ::boost::math::quantile( // Check quantile complement
61 skew_normal_distribution<RealType>(mean, scale, shape), // distribution.
62 q)), // probability complement.
63 x, // random variable.
66 skew_normal_distribution<RealType> dist (mean, scale, shape);
68 if((p < 0.999) && (q < 0.999))
69 { // We can only check this if P is not too close to 1,
70 // so that we can guarantee Q is accurate:
71 BOOST_CHECK_CLOSE_FRACTION(
72 cdf(complement(dist, x)), q, tol); // 1 - cdf
73 BOOST_CHECK_CLOSE_FRACTION(
74 quantile(dist, p), x, tol); // quantile(cdf) = x
75 BOOST_CHECK_CLOSE_FRACTION(
76 quantile(complement(dist, q)), x, tol); // quantile(complement(1 - cdf)) = x
78 } // template <class RealType>void check_skew_normal()
81 template <class RealType>
82 void test_spots(RealType)
84 // Basic sanity checks
85 RealType tolerance = 1e-4f; // 1e-4 (as %)
87 // Check some bad parameters to the distribution,
89 BOOST_CHECK_THROW(boost::math::skew_normal_distribution<RealType> nbad1(0, 0), std::domain_error); // zero sd
90 BOOST_CHECK_THROW(boost::math::skew_normal_distribution<RealType> nbad1(0, -1), std::domain_error); // negative sd
92 // Tests on extreme values of random variate x, if has numeric_limit infinity etc.
93 skew_normal_distribution<RealType> N01;
94 if(std::numeric_limits<RealType>::has_infinity)
96 BOOST_CHECK_EQUAL(pdf(N01, +std::numeric_limits<RealType>::infinity()), 0); // x = + infinity, pdf = 0
97 BOOST_CHECK_EQUAL(pdf(N01, -std::numeric_limits<RealType>::infinity()), 0); // x = - infinity, pdf = 0
98 BOOST_CHECK_EQUAL(cdf(N01, +std::numeric_limits<RealType>::infinity()), 1); // x = + infinity, cdf = 1
99 BOOST_CHECK_EQUAL(cdf(N01, -std::numeric_limits<RealType>::infinity()), 0); // x = - infinity, cdf = 0
100 BOOST_CHECK_EQUAL(cdf(complement(N01, +std::numeric_limits<RealType>::infinity())), 0); // x = + infinity, c cdf = 0
101 BOOST_CHECK_EQUAL(cdf(complement(N01, -std::numeric_limits<RealType>::infinity())), 1); // x = - infinity, c cdf = 1
102 BOOST_CHECK_THROW(boost::math::skew_normal_distribution<RealType> nbad1(std::numeric_limits<RealType>::infinity(), static_cast<RealType>(1)), std::domain_error); // +infinite mean
103 BOOST_CHECK_THROW(boost::math::skew_normal_distribution<RealType> nbad1(-std::numeric_limits<RealType>::infinity(), static_cast<RealType>(1)), std::domain_error); // -infinite mean
104 BOOST_CHECK_THROW(boost::math::skew_normal_distribution<RealType> nbad1(static_cast<RealType>(0), std::numeric_limits<RealType>::infinity()), std::domain_error); // infinite sd
107 if (std::numeric_limits<RealType>::has_quiet_NaN)
109 // No longer allow x to be NaN, then these tests should throw.
110 BOOST_CHECK_THROW(pdf(N01, +std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // x = NaN
111 BOOST_CHECK_THROW(cdf(N01, +std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // x = NaN
112 BOOST_CHECK_THROW(cdf(complement(N01, +std::numeric_limits<RealType>::quiet_NaN())), std::domain_error); // x = + infinity
113 BOOST_CHECK_THROW(quantile(N01, +std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // p = + infinity
114 BOOST_CHECK_THROW(quantile(complement(N01, +std::numeric_limits<RealType>::quiet_NaN())), std::domain_error); // p = + infinity
117 cout << "Tolerance for type " << typeid(RealType).name() << " is " << tolerance << " %" << endl;
119 // Tests where shape = 0, so same as normal tests.
120 // (These might be removed later).
122 static_cast<RealType>(5),
123 static_cast<RealType>(2),
124 static_cast<RealType>(0),
125 static_cast<RealType>(4.8),
126 static_cast<RealType>(0.46017),
127 static_cast<RealType>(1 - 0.46017),
131 static_cast<RealType>(5),
132 static_cast<RealType>(2),
133 static_cast<RealType>(0),
134 static_cast<RealType>(5.2),
135 static_cast<RealType>(1 - 0.46017),
136 static_cast<RealType>(0.46017),
140 static_cast<RealType>(5),
141 static_cast<RealType>(2),
142 static_cast<RealType>(0),
143 static_cast<RealType>(2.2),
144 static_cast<RealType>(0.08076),
145 static_cast<RealType>(1 - 0.08076),
149 static_cast<RealType>(5),
150 static_cast<RealType>(2),
151 static_cast<RealType>(0),
152 static_cast<RealType>(7.8),
153 static_cast<RealType>(1 - 0.08076),
154 static_cast<RealType>(0.08076),
158 static_cast<RealType>(-3),
159 static_cast<RealType>(5),
160 static_cast<RealType>(0),
161 static_cast<RealType>(-4.5),
162 static_cast<RealType>(0.38209),
163 static_cast<RealType>(1 - 0.38209),
167 static_cast<RealType>(-3),
168 static_cast<RealType>(5),
169 static_cast<RealType>(0),
170 static_cast<RealType>(-1.5),
171 static_cast<RealType>(1 - 0.38209),
172 static_cast<RealType>(0.38209),
176 static_cast<RealType>(-3),
177 static_cast<RealType>(5),
178 static_cast<RealType>(0),
179 static_cast<RealType>(-8.5),
180 static_cast<RealType>(0.13567),
181 static_cast<RealType>(1 - 0.13567),
185 static_cast<RealType>(-3),
186 static_cast<RealType>(5),
187 static_cast<RealType>(0),
188 static_cast<RealType>(2.5),
189 static_cast<RealType>(1 - 0.13567),
190 static_cast<RealType>(0.13567),
193 // Tests where shape != 0, specific to skew_normal distribution.
194 //void check_skew_normal(RealType mean, RealType scale, RealType shape, RealType x, RealType p, RealType q, RealType tol)
195 check_skew_normal( // 1st R example.
196 static_cast<RealType>(1.1),
197 static_cast<RealType>(2.2),
198 static_cast<RealType>(-3.3),
199 static_cast<RealType>(0.4), // x
200 static_cast<RealType>(0.733918618927874), // p == psn
201 static_cast<RealType>(1 - 0.733918618927874), // q
204 // Not sure about these yet.
205 //check_skew_normal( // 2nd R example.
206 //static_cast<RealType>(1.1),
207 //static_cast<RealType>(0.02),
208 //static_cast<RealType>(0.03),
209 //static_cast<RealType>(1.3), // x
210 //static_cast<RealType>(0.01), // p
211 //static_cast<RealType>(0.09), // q
213 //check_skew_normal( // 3nd R example.
214 //static_cast<RealType>(10.1),
215 //static_cast<RealType>(5.),
216 //static_cast<RealType>(-0.03),
217 //static_cast<RealType>(-1.3), // x
218 //static_cast<RealType>(0.01201290665838824), // p
219 //static_cast<RealType>(1. - 0.01201290665838824), // q 0.987987101
222 // Tests for PDF: we know that the normal peak value is at 1/sqrt(2*pi)
224 tolerance = boost::math::tools::epsilon<RealType>() * 5; // 5 eps as a fraction
225 BOOST_CHECK_CLOSE_FRACTION(
226 pdf(skew_normal_distribution<RealType>(), static_cast<RealType>(0)),
227 static_cast<RealType>(0.3989422804014326779399460599343818684759L), // 1/sqrt(2*pi)
229 BOOST_CHECK_CLOSE_FRACTION(
230 pdf(skew_normal_distribution<RealType>(3), static_cast<RealType>(3)),
231 static_cast<RealType>(0.3989422804014326779399460599343818684759L),
233 BOOST_CHECK_CLOSE_FRACTION(
234 pdf(skew_normal_distribution<RealType>(3, 5), static_cast<RealType>(3)),
235 static_cast<RealType>(0.3989422804014326779399460599343818684759L / 5),
239 BOOST_CHECK_CLOSE_FRACTION(
240 pdf(skew_normal_distribution<RealType>(3,5,1e-6), static_cast<RealType>(3)),
241 static_cast<RealType>(0.3989422804014326779399460599343818684759L / 5),
245 // Checks on mean, variance cumulants etc.
246 // Checks on shape ==0
248 RealType tol5 = boost::math::tools::epsilon<RealType>() * 5;
249 skew_normal_distribution<RealType> dist(8, 3);
250 RealType x = static_cast<RealType>(0.125);
252 BOOST_MATH_STD_USING // ADL of std math lib names
257 , static_cast<RealType>(8), tol5);
261 , static_cast<RealType>(9), tol5);
264 standard_deviation(dist)
265 , static_cast<RealType>(3), tol5);
269 , pdf(dist, x) / cdf(complement(dist, x)), tol5);
270 // cumulative hazard:
273 , -log(cdf(complement(dist, x))), tol5);
274 // coefficient_of_variation:
276 coefficient_of_variation(dist)
277 , standard_deviation(dist) / mean(dist), tol5);
279 BOOST_CHECK_CLOSE_FRACTION(mode(dist), static_cast<RealType>(8), 0.001);
283 , static_cast<RealType>(8), tol5);
288 , static_cast<RealType>(0), tol5);
292 , static_cast<RealType>(3), tol5);
295 kurtosis_excess(dist)
296 , static_cast<RealType>(0), tol5);
298 skew_normal_distribution<RealType> norm01(0, 1); // Test default (0, 1)
301 static_cast<RealType>(0), 0); // Mean == zero
303 skew_normal_distribution<RealType> defsd_norm01(0); // Test default (0, sd = 1)
306 static_cast<RealType>(0), 0); // Mean == zero
308 skew_normal_distribution<RealType> def_norm01; // Test default (0, sd = 1)
311 static_cast<RealType>(0), 0); // Mean == zero
314 standard_deviation(def_norm01),
315 static_cast<RealType>(1), 0); //
319 static_cast<RealType>(0), 0); // Mode == zero
322 // Skew_normal tests with shape != 0.
324 //RealType tol5 = boost::math::tools::epsilon<RealType>() * 5;
325 RealType tol100 = boost::math::tools::epsilon<RealType>() * 100;
326 RealType tol1000 = boost::math::tools::epsilon<RealType>() * 1000;
328 //skew_normal_distribution<RealType> dist(1.1, 0.02, 0.03);
330 BOOST_MATH_STD_USING // ADL of std math lib names.
332 // Test values from R = see skew_normal_drv.cpp which included the R code used.
334 skew_normal_distribution<RealType> dist(static_cast<RealType>(1.1l), static_cast<RealType>(2.2l), static_cast<RealType>(-3.3l));
336 BOOST_CHECK_CLOSE( // mean:
338 , static_cast<RealType>(-0.579908992539856825862549L), tol100);
339 BOOST_CHECK_CLOSE( // variance:
341 , static_cast<RealType>(2.0179057767837232633904L), tol100);
343 BOOST_CHECK_CLOSE( // skewness:
345 , static_cast<RealType>(-0.709854548171537509192897824663L), tol1000);
346 BOOST_CHECK_CLOSE( // kurtosis:
348 , static_cast<RealType>(3.5538752625241790601377L), tol1000);
349 BOOST_CHECK_CLOSE( // kurtosis excess:
350 kurtosis_excess(dist)
351 , static_cast<RealType>(0.5538752625241790601377L), tol1000);
354 pdf(dist, static_cast<RealType>(0.4L)),
355 static_cast<RealType>(0.294140110156599539564571L),
359 cdf(dist, static_cast<RealType>(0.4L)),
360 static_cast<RealType>(0.7339186189278737976326676452L),
364 quantile(dist, static_cast<RealType>(0.3L)),
365 static_cast<RealType>(-1.180104068086875314419247L),
371 skew_normal_distribution<RealType> dist(static_cast<RealType>(0.l), static_cast<RealType>(1.l), static_cast<RealType>(4.l));
373 cout << "pdf(dist, 0) = " << pdf(dist, 0) << ", pdf(dist, 0.45) = " << pdf(dist, 0.45) << endl;
374 // BOOST_CHECK_CLOSE(mode(dist), boost::math::constants::root_two<RealType>() / 2, tol5);
375 BOOST_CHECK_CLOSE(mode(dist), static_cast<RealType>(0.41697299497388863932L), tol1000);
381 skew_normal_distribution<RealType> dist(static_cast<RealType>(1.1l), static_cast<RealType>(0.02l), static_cast<RealType>(0.03l));
383 BOOST_CHECK_CLOSE( // mean:
385 , static_cast<RealType>(1.1004785154529557886162L), tol100);
386 BOOST_CHECK_CLOSE( // variance:
388 , static_cast<RealType>(0.00039977102296128251645L), tol100);
390 BOOST_CHECK_CLOSE( // skewness:
392 , static_cast<RealType>(5.8834811259890359782e-006L), tol1000);
393 BOOST_CHECK_CLOSE( // kurtosis:
395 , static_cast<RealType>(3.L + 9.2903475812137800239002e-008L), tol1000);
396 BOOST_CHECK_CLOSE( // kurtosis excess:
397 kurtosis_excess(dist)
398 , static_cast<RealType>(9.2903475812137800239002e-008L), tol1000);
401 skew_normal_distribution<RealType> dist(static_cast<RealType>(10.1l), static_cast<RealType>(5.l), static_cast<RealType>(-0.03l));
402 BOOST_CHECK_CLOSE( // mean:
404 , static_cast<RealType>(9.9803711367610528459485937L), tol100);
405 BOOST_CHECK_CLOSE( // variance:
407 , static_cast<RealType>(24.98568893508015727823L), tol100);
409 BOOST_CHECK_CLOSE( // skewness:
411 , static_cast<RealType>(-5.8834811259890359782085e-006L), tol1000);
412 BOOST_CHECK_CLOSE( // kurtosis:
414 , static_cast<RealType>(3.L + 9.2903475812137800239002e-008L), tol1000);
415 BOOST_CHECK_CLOSE( // kurtosis excess:
416 kurtosis_excess(dist)
417 , static_cast<RealType>(9.2903475812137800239002e-008L), tol1000);
420 skew_normal_distribution<RealType> dist(static_cast<RealType>(-10.1l), static_cast<RealType>(5.l), static_cast<RealType>(30.l));
421 BOOST_CHECK_CLOSE( // mean:
423 , static_cast<RealType>(-6.11279169674138408531365L), 2 * tol100);
424 BOOST_CHECK_CLOSE( // variance:
426 , static_cast<RealType>(9.10216994642554914628242L), tol100 * 2);
428 BOOST_CHECK_CLOSE( // skewness:
430 , static_cast<RealType>(0.99072425443686904424L), tol1000);
431 BOOST_CHECK_CLOSE( // kurtosis:
433 , static_cast<RealType>(3.L + 0.8638862008406084244563L), tol1000);
434 BOOST_CHECK_CLOSE( // kurtosis excess:
435 kurtosis_excess(dist)
436 , static_cast<RealType>(0.8638862008406084244563L), tol1000);
443 } // template <class RealType>void test_spots(RealType)
445 int test_main(int, char* [])
449 using boost::math::skew_normal;
450 using boost::math::skew_normal_distribution;
452 //int precision = 17; // std::numeric_limits<double::max_digits10;
453 double tolfeweps = numeric_limits<double>::epsilon() * 5;
454 //double tol6decdigits = numeric_limits<float>::epsilon() * 2;
455 // Check that can generate skew_normal distribution using the two convenience methods:
456 boost::math::skew_normal w12(1., 2); // Using typedef.
457 boost::math::skew_normal_distribution<> w01; // Use default unity values for mean and scale.
458 // Note NOT myn01() as the compiler will interpret as a function!
460 // Checks on constructors.
461 // Default parameters.
462 BOOST_CHECK_EQUAL(w01.location(), 0);
463 BOOST_CHECK_EQUAL(w01.scale(), 1);
464 BOOST_CHECK_EQUAL(w01.shape(), 0);
466 skew_normal_distribution<> w23(2., 3); // Using default RealType double.
467 BOOST_CHECK_EQUAL(w23.scale(), 3);
468 BOOST_CHECK_EQUAL(w23.shape(), 0);
470 skew_normal_distribution<> w123(1., 2., 3.); // Using default RealType double.
471 BOOST_CHECK_EQUAL(w123.location(), 1.);
472 BOOST_CHECK_EQUAL(w123.scale(), 2.);
473 BOOST_CHECK_EQUAL(w123.shape(), 3.);
475 BOOST_CHECK_CLOSE_FRACTION(mean(w01), static_cast<double>(0), tolfeweps); // Default mean == zero
476 BOOST_CHECK_CLOSE_FRACTION(scale(w01), static_cast<double>(1), tolfeweps); // Default scale == unity
478 // Basic sanity-check spot values for all floating-point types..
479 // (Parameter value, arbitrarily zero, only communicates the floating point type).
480 test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 %
481 test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 %
482 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
483 test_spots(0.0L); // Test long double.
484 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
485 test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
488 std::cout << "<note>The long double tests have been disabled on this platform "
489 "either because the long double overloads of the usual math functions are "
490 "not available at all, or because they are too inaccurate for these tests "
491 "to pass.</note>" << std::cout;
495 } // int test_main(int, char* [])