Subject: [Boost-commit] svn:boost r76718 - in sandbox/math/libs: . math math/doc math/example math/test
From: pbristow_at_[hidden]
Date: 2012-01-27 12:08:57
Author: pbristow
Date: 2012-01-27 12:08:55 EST (Fri, 27 Jan 2012)
New Revision: 76718
URL: http://svn.boost.org/trac/boost/changeset/76718
Log:
1st commit of skew_normal to boost sandbox, but failures in mode TODO.
Added:
sandbox/math/libs/
sandbox/math/libs/math/
sandbox/math/libs/math/doc/
sandbox/math/libs/math/example/
sandbox/math/libs/math/example/owens_t_drv.cpp (contents, props changed)
sandbox/math/libs/math/example/skew_normal_drv.cpp (contents, props changed)
sandbox/math/libs/math/test/
sandbox/math/libs/math/test/test_owens_t.cpp (contents, props changed)
sandbox/math/libs/math/test/test_skew_normal.cpp (contents, props changed)
Added: sandbox/math/libs/math/example/owens_t_drv.cpp
==============================================================================
--- (empty file)
+++ sandbox/math/libs/math/example/owens_t_drv.cpp 2012-01-27 12:08:55 EST (Fri, 27 Jan 2012)
@@ -0,0 +1,104 @@
+
+#include "owens_t.hpp"
+#include <iostream>
+
+int main()
+{
+ double h = 0.0,a;
+ std::cout << std::setprecision(20);
+
+ static const double a_vec[] = {
+ 0.5000000000000000E+00,
+ 0.1000000000000000E+01,
+ 0.2000000000000000E+01,
+ 0.3000000000000000E+01,
+ 0.5000000000000000E+00,
+ 0.1000000000000000E+01,
+ 0.2000000000000000E+01,
+ 0.3000000000000000E+01,
+ 0.5000000000000000E+00,
+ 0.1000000000000000E+01,
+ 0.2000000000000000E+01,
+ 0.3000000000000000E+01,
+ 0.5000000000000000E+00,
+ 0.1000000000000000E+01,
+ 0.2000000000000000E+01,
+ 0.3000000000000000E+01,
+ 0.5000000000000000E+00,
+ 0.1000000000000000E+01,
+ 0.2000000000000000E+01,
+ 0.3000000000000000E+01,
+ 0.1000000000000000E+02,
+ 0.1000000000000000E+03 };
+
+ static const double h_vec[] = {
+ 0.1000000000000000E+01,
+ 0.1000000000000000E+01,
+ 0.1000000000000000E+01,
+ 0.1000000000000000E+01,
+ 0.5000000000000000E+00,
+ 0.5000000000000000E+00,
+ 0.5000000000000000E+00,
+ 0.5000000000000000E+00,
+ 0.2500000000000000E+00,
+ 0.2500000000000000E+00,
+ 0.2500000000000000E+00,
+ 0.2500000000000000E+00,
+ 0.1250000000000000E+00,
+ 0.1250000000000000E+00,
+ 0.1250000000000000E+00,
+ 0.1250000000000000E+00,
+ 0.7812500000000000E-02,
+ 0.7812500000000000E-02,
+ 0.7812500000000000E-02,
+ 0.7812500000000000E-02,
+ 0.7812500000000000E-02,
+ 0.7812500000000000E-02 };
+
+ static const double t_vec[] = {
+ 0.4306469112078537E-01,
+ 0.6674188216570097E-01,
+ 0.7846818699308410E-01,
+ 0.7929950474887259E-01,
+ 0.6448860284750376E-01,
+ 0.1066710629614485E+00,
+ 0.1415806036539784E+00,
+ 0.1510840430760184E+00,
+ 0.7134663382271778E-01,
+ 0.1201285306350883E+00,
+ 0.1666128410939293E+00,
+ 0.1847501847929859E+00,
+ 0.7317273327500385E-01,
+ 0.1237630544953746E+00,
+ 0.1737438887583106E+00,
+ 0.1951190307092811E+00,
+ 0.7378938035365546E-01,
+ 0.1249951430754052E+00,
+ 0.1761984774738108E+00,
+ 0.1987772386442824E+00,
+ 0.2340886964802671E+00,
+ 0.2479460829231492E+00 };
+
+ for(unsigned i = 0; i != 22; ++i)
+ {
+ h = h_vec[i];
+ a = a_vec[i];
+ const double t = boost::math::owens_t(h, a);
+ std::cout << "h=" << h << "\ta=" << a << "\tcomp="
+ << t << "\ttab=" << t_vec[i]
+ << "\tdiff=" << fabs(t_vec[i]-t) << std::endl;;
+ }
+
+ return 0;
+}
+
+
+// EOF
+
+
+
+
+
+
+
+
Added: sandbox/math/libs/math/example/skew_normal_drv.cpp
==============================================================================
--- (empty file)
+++ sandbox/math/libs/math/example/skew_normal_drv.cpp 2012-01-27 12:08:55 EST (Fri, 27 Jan 2012)
@@ -0,0 +1,213 @@
+// (C) Benjamin Sobotta 2012
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+
+#ifdef _MSC_VER
+# pragma warning (disable : 4512) // assignment operator could not be generated
+# pragma warning (disable : 4127) // conditional expression is constant.
+#endif
+
+#include <boost/math/distributions/skew_normal.hpp>
+using boost::math::skew_normal_distribution;
+using boost::math::skew_normal;
+#include <iostream>
+#include <cmath>
+#include <utility>
+
+void check(const double loc, const double sc, const double sh,
+ const double * const cumulants, const std::pair<double, double> qu,
+ const double x, const double tpdf, const double tcdf)
+{
+ using namespace boost::math;
+
+ skew_normal D(loc, sc, sh);
+
+ const double sk = cumulants[2] / (std::pow(cumulants[1], 1.5));
+ const double kt = cumulants[3] / (cumulants[1] * cumulants[1]);
+
+ // checks against tabulated values
+ std::cout << "mean: table=" << cumulants[0] << "\tcompute=" << mean(D) << "\tdiff=" << fabs(cumulants[0]-mean(D)) << std::endl;
+ std::cout << "var: table=" << cumulants[1] << "\tcompute=" << variance(D) << "\tdiff=" << fabs(cumulants[0]-variance(D)) << std::endl;
+ std::cout << "skew: table=" << sk << "\tcompute=" << skewness(D) << "\tdiff=" << fabs(sk-skewness(D)) << std::endl;
+ std::cout << "kur.: table=" << kt << "\tcompute=" << kurtosis_excess(D) << "\tdiff=" << fabs(kt-kurtosis_excess(D)) << std::endl;
+
+ const double q = quantile(D, qu.first);
+ const double cq = quantile(complement(D, qu.first));
+
+ std::cout << "quantile(" << qu.first << "): table=" << qu.second << "\tcompute=" << q << "\tdiff=" << fabs(qu.second-q) << std::endl;
+
+ // consistency
+ std::cout << "cdf(quantile)=" << cdf(D, q) << "\tp=" << qu.first << "\tdiff=" << fabs(qu.first-cdf(D, q)) << std::endl;
+ std::cout << "ccdf(cquantile)=" << cdf(complement(D,cq)) << "\tp=" << qu.first << "\tdiff=" << fabs(qu.first-cdf(complement(D,cq))) << std::endl;
+
+ // PDF & CDF
+ std::cout << "pdf(" << x << "): table=" << tpdf << "\tcompute=" << pdf(D,x) << "\tdiff=" << fabs(tpdf-pdf(D,x)) << std::endl;
+ std::cout << "cdf(" << x << "): table=" << tcdf << "\tcompute=" << cdf(D,x) << "\tdiff=" << fabs(tcdf-cdf(D,x)) << std::endl;
+ std::cout << "================================\n";
+}
+
+int main()
+{
+ using namespace boost::math;
+
+ double sc = 0.0,loc,sh,x,dsn,qsn,psn,p;
+ std::cout << std::setprecision(20);
+
+ double cumulants[4];
+
+
+ /* R:
+ > install.packages("sn")
+ Warning in install.packages("sn") :
+ 'lib = "/usr/lib64/R/library"' is not writable
+ Would you like to create a personal library
+ '~/R/x86_64-unknown-linux-gnu-library/2.12'
+ to install packages into? (y/n) y
+ --- Please select a CRAN mirror for use in this session ---
+ Loading Tcl/Tk interface ... done
+ also installing the dependency ‘mnormt’
+
+ trying URL 'http://mirrors.dotsrc.org/cran/src/contrib/mnormt_1.4-5.tar.gz'
+ Content type 'application/x-gzip' length 34049 bytes (33 Kb)
+ opened URL
+ ==================================================
+ downloaded 33 Kb
+
+ trying URL 'http://mirrors.dotsrc.org/cran/src/contrib/sn_0.4-17.tar.gz'
+ Content type 'application/x-gzip' length 65451 bytes (63 Kb)
+ opened URL
+ ==================================================
+ downloaded 63 Kb
+
+
+ > library(sn)
+ > options(digits=22)
+
+
+ > sn.cumulants(1.1, 2.2, -3.3)
+ [1] -0.5799089925398568 2.0179057767837230 -2.0347951542374196
+ [4] 2.2553488991015072
+ > qsn(0.3, 1.1, 2.2, -3.3)
+ [1] -1.180104068086876
+ > psn(0.4, 1.1, 2.2, -3.3)
+ [1] 0.733918618927874
+ > dsn(0.4, 1.1, 2.2, -3.3)
+ [1] 0.2941401101565995
+
+ */
+
+ //1 st
+ loc = 1.1; sc = 2.2; sh = -3.3;
+ cumulants[0] = -0.5799089925398568;
+ cumulants[1] = 2.0179057767837230;
+ cumulants[2] = -2.0347951542374196;
+ cumulants[3] = 2.2553488991015072;
+ x = 0.4;
+ p = 0.3;
+ qsn = -1.180104068086876;
+ psn = 0.733918618927874;
+ dsn = 0.2941401101565995;
+
+ check(loc, sc, sh, cumulants, std::make_pair(p,qsn), x, dsn, psn);
+
+ /* R:
+
+ > sn.cumulants(1.1, .02, .03)
+ [1] 1.1004785154529559e+00 3.9977102296128255e-04 4.7027439329779991e-11
+ [4] 1.4847542790693825e-14
+ > qsn(0.01, 1.1, .02, .03)
+ [1] 1.053964962950150
+ > psn(1.3, 1.1, .02, .03)
+ [1] 1
+ > dsn(1.3, 1.1, .02, .03)
+ [1] 4.754580380601393e-21
+
+ */
+
+ // 2nd
+ loc = 1.1; sc = .02; sh = .03;
+ cumulants[0] = 1.1004785154529559e+00;
+ cumulants[1] = 3.9977102296128255e-04;
+ cumulants[2] = 4.7027439329779991e-11;
+ cumulants[3] = 1.4847542790693825e-14;
+ x = 1.3;
+ p = 0.01;
+ qsn = 1.053964962950150;
+ psn = 1;
+ dsn = 4.754580380601393e-21;
+
+ check(loc, sc, sh, cumulants, std::make_pair(p,qsn), x, dsn, psn);
+
+ /* R:
+
+ > sn.cumulants(10.1, 5, -.03)
+ [1] 9.980371136761052e+00 2.498568893508016e+01 -7.348037395278123e-04
+ [4] 5.799821402614775e-05
+ > qsn(.8, 10.1, 5, -.03)
+ [1] 14.18727407491953
+ > psn(-1.3, 10.1, 5, -.03)
+ [1] 0.01201290665838824
+ > dsn(-1.3, 10.1, 5, -.03)
+ [1] 0.006254346348897927
+
+
+ */
+
+ // 3rd
+ loc = 10.1; sc = 5; sh = -.03;
+ cumulants[0] = 9.980371136761052e+00;
+ cumulants[1] = 2.498568893508016e+01;
+ cumulants[2] = -7.348037395278123e-04;
+ cumulants[3] = 5.799821402614775e-05;
+ x = -1.3;
+ p = 0.8;
+ qsn = 14.18727407491953;
+ psn = 0.01201290665838824;
+ dsn = 0.006254346348897927;
+
+ check(loc, sc, sh, cumulants, std::make_pair(p,qsn), x, dsn, psn);
+
+
+ /* R:
+
+ > sn.cumulants(-10.1, 5, 30)
+ [1] -6.112791696741384 9.102169946425548 27.206345266148194 71.572537838642916
+ > qsn(.8, -10.1, 5, 30)
+ [1] -3.692242172277
+ > psn(-1.3, -10.1, 5, 30)
+ [1] 0.921592193425035
+ > dsn(-1.3, -10.1, 5, 30)
+ [1] 0.0339105445232089
+
+ */
+
+ // 4th
+ loc = -10.1; sc = 5; sh = 30;
+ cumulants[0] = -6.112791696741384;
+ cumulants[1] = 9.102169946425548;
+ cumulants[2] = 27.206345266148194;
+ cumulants[3] = 71.572537838642916;
+ x = -1.3;
+ p = 0.8;
+ qsn = -3.692242172277;
+ psn = 0.921592193425035;
+ dsn = 0.0339105445232089;
+
+ check(loc, sc, sh, cumulants, std::make_pair(p,qsn), x, dsn, psn);
+
+ return 0;
+}
+
+
+// EOF
+
+
+
+
+
+
+
+
Added: sandbox/math/libs/math/test/test_owens_t.cpp
==============================================================================
--- (empty file)
+++ sandbox/math/libs/math/test/test_owens_t.cpp 2012-01-27 12:08:55 EST (Fri, 27 Jan 2012)
@@ -0,0 +1,153 @@
+// test_owens_t.cpp
+
+// Copyright Paul A. Bristow 2012.
+// Copyright Benjamin Sobotta 2012.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Tested using some 30 decimal digit accuracy values from:
+// Fast and accurate calculation of Owen's T-function
+// Mike Patefield, and David Tandy
+// Journal of Statistical Software, 5 (5), 1-25 (2000).
+// http://www.jstatsoft.org/v05/a05/paper Table 3, page 15
+// Values of T(h,a) accurate to thirty figures were calculated using 128 bit arithmetic by
+// evaluating (9) with m = 48, the summation over k being continued until additional terms did
+// not alter the result. The resultant values Tacc(h,a) say, were validated by evaluating (8) with
+// m = 48 (i.e. 96 point Gaussian quadrature).
+
+#ifdef _MSC_VER
+# pragma warning (disable : 4127) // conditional expression is constant
+# pragma warning (disable : 4305) // 'initializing' : truncation from 'double' to 'const float'
+// ?? TODO get rid of these warnings?
+#endif
+
+#include <boost/math/concepts/real_concept.hpp> // for real_concept.
+using ::boost::math::concepts::real_concept;
+
+#include <boost/math/special_functions/owens_t.hpp> // for owens_t function.
+using boost::math::owens_t;
+
+#include <boost/test/test_exec_monitor.hpp> // for test_main.
+#include <boost/test/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE and BOOST_CHECK_CLOSE_fraction.
+
+#include <iostream>
+using std::cout;
+using std::endl;
+#include <limits>
+using std::numeric_limits;
+
+template <class RealType>
+void test_spot(
+ RealType h, //
+ RealType a, //
+ RealType tol) // Test tolerance
+{
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(h, a), 3.89119302347013668966224771378e-2L, tol);
+}
+
+
+template <class RealType> // Any floating-point type RealType.
+void test_spots(RealType)
+{
+ // Basic sanity checks, test data is as accurate as long double,
+ // so set tolerance to a few epsilon expressed as a fraction.
+ RealType tolerance = boost::math::tools::epsilon<RealType>() * 30; // most OK with 3 eps tolerance.
+ cout << "Tolerance = " << tolerance << "." << endl;
+
+ using ::boost::math::owens_t;
+ using namespace std; // ADL of std names.
+
+ // Checks of six sub-methods T1 to T6.
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.0625L), static_cast<RealType>(0.25L)), static_cast<RealType>(3.89119302347013668966224771378e-2L), tolerance); // T1
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(6.5L), static_cast<RealType>(0.4375L)), static_cast<RealType>(2.00057730485083154100907167685E-11L), tolerance); // T2
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(7L), static_cast<RealType>( 0.96875L)), static_cast<RealType>(6.39906271938986853083219914429E-13L), tolerance); // T3
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(4.78125L), static_cast<RealType>(0.0625L)), static_cast<RealType>(1.06329748046874638058307112826E-7L), tolerance); // T4
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(2.L), static_cast<RealType>(0.5L)), static_cast<RealType>(8.62507798552150713113488319155E-3L), tolerance); // T5
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(1.L), static_cast<RealType>(0.9999975L)), static_cast<RealType>(6.67418089782285927715589822405E-2L), tolerance); // T6
+ //BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(L), static_cast<RealType>(L)), static_cast<RealType>(L), tolerance);
+
+ // Mathematica spot values.
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.4375L), static_cast<RealType>(6.5L)), static_cast<RealType>(0.1654013012544939L), tolerance);
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.96875L), static_cast<RealType>(7.L)), static_cast<RealType>(0.0831674847460297L), tolerance);
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.0625L), static_cast<RealType>(4.78125L)), static_cast<RealType>(0.2157118581989799L), tolerance);
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.5L), static_cast<RealType>(2.L)), static_cast<RealType>(0.1415806036539784L), tolerance);
+
+
+// BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(L), static_cast<RealType>(L)), static_cast<RealType>(L), tolerance);
+
+ // Spots values using Mathematica but only 17 decimal digits.
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(6.5L), static_cast<RealType>(0.4375L)), static_cast<RealType>(2.000577304850829E-11L), tolerance);
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.4375L), static_cast<RealType>(6.5L)), static_cast<RealType>(0.1654013012544939L), tolerance);
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(7.L), static_cast<RealType>(0.96875L)), static_cast<RealType>(6.399062719389912e-13L), tolerance);
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.96875L), static_cast<RealType>(7.L)), static_cast<RealType>(0.0831674847460297L), tolerance);
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(4.78125L), static_cast<RealType>(0.0625L)), static_cast<RealType>(1.063297480468747e-7L), tolerance);
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.0625L), static_cast<RealType>(4.78125L)), static_cast<RealType>(0.2157118581989799L), tolerance);
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(2.L), static_cast<RealType>(0.5L)), static_cast<RealType>(0.00862507798552151L), tolerance);
+ BOOST_CHECK_CLOSE_FRACTION(owens_t(static_cast<RealType>(0.5L), static_cast<RealType>(2L)), static_cast<RealType>(0.1415806036539784L), tolerance);
+
+
+} // template <class RealType>void test_spots(RealType)
+
+int test_main(int, char* [])
+{
+ BOOST_MATH_CONTROL_FP;
+
+ // Basic sanity-check spot values.
+
+ // (Parameter value, arbitrarily zero, only communicates the floating point type).
+ test_spots(0.0F); // Test float.
+ test_spots(0.0); // Test double.
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+ test_spots(0.0L); // Test long double.
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
+ test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
+#endif
+#endif
+ return 0;
+} // int test_main(int, char* [])
+
+/*
+
+Output:
+------ Rebuild All started: Project: test_owens_t_sandbox, Configuration: Release Win32 ------
+ test_owens_t.cpp
+ Generating code
+ Finished generating code
+ test_owens_t_sandbox.vcxproj -> J:\Cpp\MathToolkit\test\Math_test\Release\test_owens_t_sandbox.exe
+ Running 1 test case...
+ Platform: Win32
+ Compiler: Microsoft Visual C++ version 10.0
+ STL : Dinkumware standard library version 520
+ Boost : 1.49.0
+ Tolerance = 3.57628e-007.
+ Tolerance = 6.66134e-016.
+ Tolerance = 6.66134e-016.
+ Tolerance = 6.66134e-016.
+
+ *** No errors detected
+========== Rebuild All: 1 succeeded, 0 failed, 0 skipped ==========
+
+
+------ Build started: Project: test_owens_t_sandbox, Configuration: Debug Win32 ------
+ test_owens_t.cpp
+ test_owens_t_sandbox.vcxproj -> J:\Cpp\MathToolkit\test\Math_test\Debug\test_owens_t_sandbox.exe
+ Running 1 test case...
+ Platform: Win32
+ Compiler: Microsoft Visual C++ version 10.0
+ STL : Dinkumware standard library version 520
+ Boost : 1.49.0
+ Tolerance = 3.57628e-006.
+ Tolerance = 6.66134e-015.
+ Tolerance = 6.66134e-015.
+ Tolerance = 6.66134e-015.
+
+ *** No errors detected
+========== Build: 1 succeeded, 0 failed, 0 up-to-date, 0 skipped ==========
+
+*/
+
+
+
Added: sandbox/math/libs/math/test/test_skew_normal.cpp
==============================================================================
--- (empty file)
+++ sandbox/math/libs/math/test/test_skew_normal.cpp 2012-01-27 12:08:55 EST (Fri, 27 Jan 2012)
@@ -0,0 +1,504 @@
+// Copyright Paul A. Bristow 2012.
+// Copyright John Maddock 2012.
+// Copyright Benjamin Sobotta 2012
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifdef _MSC_VER
+# pragma warning (disable : 4127) // conditional expression is constant.
+# pragma warning (disable : 4305) // 'initializing' : truncation from 'double' to 'const float'.
+# pragma warning (disable : 4310) // cast truncates constant value.
+# pragma warning (disable : 4512) // assignment operator could not be generated.
+#endif
+
+//#include <pch.hpp> // include directory libs/math/src/tr1/ is needed.
+
+#include <boost/math/concepts/real_concept.hpp> // for real_concept
+#include <boost/test/test_exec_monitor.hpp> // Boost.Test
+#include <boost/test/floating_point_comparison.hpp>
+
+#include <boost/math/distributions/skew_normal.hpp>
+using boost::math::skew_normal_distribution;
+using boost::math::skew_normal;
+
+#include <iostream>
+using std::cout;
+using std::endl;
+using std::setprecision;
+#include <limits>
+using std::numeric_limits;
+
+template <class RealType>
+void check_skew_normal(RealType mean, RealType scale, RealType shape, RealType x, RealType p, RealType q, RealType tol)
+{
+ using boost::math::skew_normal_distribution;
+
+ BOOST_CHECK_CLOSE_FRACTION(
+ ::boost::math::cdf( // Check cdf
+ skew_normal_distribution<RealType>(mean, scale, shape), // distribution.
+ x), // random variable.
+ p, // probability.
+ tol); // tolerance.
+ BOOST_CHECK_CLOSE_FRACTION(
+ ::boost::math::cdf( // Check cdf complement
+ complement(
+ skew_normal_distribution<RealType>(mean, scale, shape), // distribution.
+ x)), // random variable.
+ q, // probability complement.
+ tol); // %tolerance.
+ BOOST_CHECK_CLOSE_FRACTION(
+ ::boost::math::quantile( // Check quantile
+ skew_normal_distribution<RealType>(mean, scale, shape), // distribution.
+ p), // probability.
+ x, // random variable.
+ tol); // tolerance.
+ BOOST_CHECK_CLOSE_FRACTION(
+ ::boost::math::quantile( // Check quantile complement
+ complement(
+ skew_normal_distribution<RealType>(mean, scale, shape), // distribution.
+ q)), // probability complement.
+ x, // random variable.
+ tol); // tolerance.
+
+ skew_normal_distribution<RealType> dist (mean, scale, shape);
+
+ if((p < 0.999) && (q < 0.999))
+ { // We can only check this if P is not too close to 1,
+ // so that we can guarantee Q is accurate:
+ BOOST_CHECK_CLOSE_FRACTION(
+ cdf(complement(dist, x)), q, tol); // 1 - cdf
+ BOOST_CHECK_CLOSE_FRACTION(
+ quantile(dist, p), x, tol); // quantile(cdf) = x
+ BOOST_CHECK_CLOSE_FRACTION(
+ quantile(complement(dist, q)), x, tol); // quantile(complement(1 - cdf)) = x
+ }
+} // template <class RealType>void check_skew_normal()
+
+
+template <class RealType>
+void test_spots(RealType)
+{
+ // Basic sanity checks
+ RealType tolerance = 1e-4f; // 1e-4 (as %)
+
+ // Check some bad parameters to the distribution,
+
+ BOOST_CHECK_THROW(boost::math::skew_normal_distribution<RealType> nbad1(0, 0), std::domain_error); // zero sd
+ BOOST_CHECK_THROW(boost::math::skew_normal_distribution<RealType> nbad1(0, -1), std::domain_error); // negative sd
+
+ // Tests on extreme values of random variate x, if has numeric_limit infinity etc.
+ skew_normal_distribution<RealType> N01;
+ if(std::numeric_limits<RealType>::has_infinity)
+ {
+ BOOST_CHECK_EQUAL(pdf(N01, +std::numeric_limits<RealType>::infinity()), 0); // x = + infinity, pdf = 0
+ BOOST_CHECK_EQUAL(pdf(N01, -std::numeric_limits<RealType>::infinity()), 0); // x = - infinity, pdf = 0
+ BOOST_CHECK_EQUAL(cdf(N01, +std::numeric_limits<RealType>::infinity()), 1); // x = + infinity, cdf = 1
+ BOOST_CHECK_EQUAL(cdf(N01, -std::numeric_limits<RealType>::infinity()), 0); // x = - infinity, cdf = 0
+ BOOST_CHECK_EQUAL(cdf(complement(N01, +std::numeric_limits<RealType>::infinity())), 0); // x = + infinity, c cdf = 0
+ BOOST_CHECK_EQUAL(cdf(complement(N01, -std::numeric_limits<RealType>::infinity())), 1); // x = - infinity, c cdf = 1
+ BOOST_CHECK_THROW(boost::math::skew_normal_distribution<RealType> nbad1(std::numeric_limits<RealType>::infinity(), static_cast<RealType>(1)), std::domain_error); // +infinite mean
+ BOOST_CHECK_THROW(boost::math::skew_normal_distribution<RealType> nbad1(-std::numeric_limits<RealType>::infinity(), static_cast<RealType>(1)), std::domain_error); // -infinite mean
+ BOOST_CHECK_THROW(boost::math::skew_normal_distribution<RealType> nbad1(static_cast<RealType>(0), std::numeric_limits<RealType>::infinity()), std::domain_error); // infinite sd
+ }
+
+ if (std::numeric_limits<RealType>::has_quiet_NaN)
+ {
+ // No longer allow x to be NaN, then these tests should throw.
+ BOOST_CHECK_THROW(pdf(N01, +std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // x = NaN
+ BOOST_CHECK_THROW(cdf(N01, +std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // x = NaN
+ BOOST_CHECK_THROW(cdf(complement(N01, +std::numeric_limits<RealType>::quiet_NaN())), std::domain_error); // x = + infinity
+ BOOST_CHECK_THROW(quantile(N01, +std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // p = + infinity
+ BOOST_CHECK_THROW(quantile(complement(N01, +std::numeric_limits<RealType>::quiet_NaN())), std::domain_error); // p = + infinity
+ }
+
+ cout << "Tolerance for type " << typeid(RealType).name() << " is " << tolerance << " %" << endl;
+
+ // Tests where shape = 0, so same as normal tests.
+ // (These might be removed later).
+ check_skew_normal(
+ static_cast<RealType>(5),
+ static_cast<RealType>(2),
+ static_cast<RealType>(0),
+ static_cast<RealType>(4.8),
+ static_cast<RealType>(0.46017),
+ static_cast<RealType>(1 - 0.46017),
+ tolerance);
+
+ check_skew_normal(
+ static_cast<RealType>(5),
+ static_cast<RealType>(2),
+ static_cast<RealType>(0),
+ static_cast<RealType>(5.2),
+ static_cast<RealType>(1 - 0.46017),
+ static_cast<RealType>(0.46017),
+ tolerance);
+
+ check_skew_normal(
+ static_cast<RealType>(5),
+ static_cast<RealType>(2),
+ static_cast<RealType>(0),
+ static_cast<RealType>(2.2),
+ static_cast<RealType>(0.08076),
+ static_cast<RealType>(1 - 0.08076),
+ tolerance);
+
+ check_skew_normal(
+ static_cast<RealType>(5),
+ static_cast<RealType>(2),
+ static_cast<RealType>(0),
+ static_cast<RealType>(7.8),
+ static_cast<RealType>(1 - 0.08076),
+ static_cast<RealType>(0.08076),
+ tolerance);
+
+ check_skew_normal(
+ static_cast<RealType>(-3),
+ static_cast<RealType>(5),
+ static_cast<RealType>(0),
+ static_cast<RealType>(-4.5),
+ static_cast<RealType>(0.38209),
+ static_cast<RealType>(1 - 0.38209),
+ tolerance);
+
+ check_skew_normal(
+ static_cast<RealType>(-3),
+ static_cast<RealType>(5),
+ static_cast<RealType>(0),
+ static_cast<RealType>(-1.5),
+ static_cast<RealType>(1 - 0.38209),
+ static_cast<RealType>(0.38209),
+ tolerance);
+
+ check_skew_normal(
+ static_cast<RealType>(-3),
+ static_cast<RealType>(5),
+ static_cast<RealType>(0),
+ static_cast<RealType>(-8.5),
+ static_cast<RealType>(0.13567),
+ static_cast<RealType>(1 - 0.13567),
+ tolerance);
+
+ check_skew_normal(
+ static_cast<RealType>(-3),
+ static_cast<RealType>(5),
+ static_cast<RealType>(0),
+ static_cast<RealType>(2.5),
+ static_cast<RealType>(1 - 0.13567),
+ static_cast<RealType>(0.13567),
+ tolerance);
+
+ // Tests where shape != 0, specific to skew_normal distribution.
+ //void check_skew_normal(RealType mean, RealType scale, RealType shape, RealType x, RealType p, RealType q, RealType tol)
+ check_skew_normal( // 1st R example.
+ static_cast<RealType>(1.1),
+ static_cast<RealType>(2.2),
+ static_cast<RealType>(-3.3),
+ static_cast<RealType>(0.4), // x
+ static_cast<RealType>(0.733918618927874), // p == psn
+ static_cast<RealType>(1 - 0.733918618927874), // q
+ tolerance);
+
+ // Not sure about these yet.
+ //check_skew_normal( // 2nd R example.
+ //static_cast<RealType>(1.1),
+ //static_cast<RealType>(0.02),
+ //static_cast<RealType>(0.03),
+ //static_cast<RealType>(1.3), // x
+ //static_cast<RealType>(0.01), // p
+ //static_cast<RealType>(0.09), // q
+ //tolerance);
+ //check_skew_normal( // 3nd R example.
+ //static_cast<RealType>(10.1),
+ //static_cast<RealType>(5.),
+ //static_cast<RealType>(-0.03),
+ //static_cast<RealType>(-1.3), // x
+ //static_cast<RealType>(0.01201290665838824), // p
+ //static_cast<RealType>(1. - 0.01201290665838824), // q 0.987987101
+ //tolerance);
+
+ // Tests for PDF: we know that the normal peak value is at 1/sqrt(2*pi)
+ //
+ tolerance = boost::math::tools::epsilon<RealType>() * 5; // 5 eps as a fraction
+ BOOST_CHECK_CLOSE_FRACTION(
+ pdf(skew_normal_distribution<RealType>(), static_cast<RealType>(0)),
+ static_cast<RealType>(0.3989422804014326779399460599343818684759L), // 1/sqrt(2*pi)
+ tolerance);
+ BOOST_CHECK_CLOSE_FRACTION(
+ pdf(skew_normal_distribution<RealType>(3), static_cast<RealType>(3)),
+ static_cast<RealType>(0.3989422804014326779399460599343818684759L),
+ tolerance);
+ BOOST_CHECK_CLOSE_FRACTION(
+ pdf(skew_normal_distribution<RealType>(3, 5), static_cast<RealType>(3)),
+ static_cast<RealType>(0.3989422804014326779399460599343818684759L / 5),
+ tolerance);
+
+ // Shape != 0.
+ BOOST_CHECK_CLOSE_FRACTION(
+ pdf(skew_normal_distribution<RealType>(1.1, 0.02, 0.03), static_cast<RealType>(3)),
+ static_cast<RealType>(0),
+ tolerance);
+
+
+ // Checks on mean, variance cumulants etc.
+ // Checks on shape ==0
+
+ RealType tol5 = boost::math::tools::epsilon<RealType>() * 5;
+ skew_normal_distribution<RealType> dist(8, 3);
+ RealType x = static_cast<RealType>(0.125);
+
+ BOOST_MATH_STD_USING // ADL of std math lib names
+
+ // mean:
+ BOOST_CHECK_CLOSE(
+ mean(dist)
+ , static_cast<RealType>(8), tol5);
+ // variance:
+ BOOST_CHECK_CLOSE(
+ variance(dist)
+ , static_cast<RealType>(9), tol5);
+ // std deviation:
+ BOOST_CHECK_CLOSE(
+ standard_deviation(dist)
+ , static_cast<RealType>(3), tol5);
+ // hazard:
+ BOOST_CHECK_CLOSE(
+ hazard(dist, x)
+ , pdf(dist, x) / cdf(complement(dist, x)), tol5);
+ // cumulative hazard:
+ BOOST_CHECK_CLOSE(
+ chf(dist, x)
+ , -log(cdf(complement(dist, x))), tol5);
+ // coefficient_of_variation:
+ BOOST_CHECK_CLOSE(
+ coefficient_of_variation(dist)
+ , standard_deviation(dist) / mean(dist), tol5);
+ // mode:
+ BOOST_CHECK_CLOSE_FRACTION(mode(dist), static_cast<RealType>(8), 0.001);
+
+ BOOST_CHECK_CLOSE(
+ median(dist)
+ , static_cast<RealType>(8), tol5);
+
+ // skewness:
+ BOOST_CHECK_CLOSE(
+ skewness(dist)
+ , static_cast<RealType>(0), tol5);
+ // kurtosis:
+ BOOST_CHECK_CLOSE(
+ kurtosis(dist)
+ , static_cast<RealType>(3), tol5);
+ // kurtosis excess:
+ BOOST_CHECK_CLOSE(
+ kurtosis_excess(dist)
+ , static_cast<RealType>(0), tol5);
+
+ skew_normal_distribution<RealType> norm01(0, 1); // Test default (0, 1)
+ BOOST_CHECK_CLOSE(
+ mean(norm01),
+ static_cast<RealType>(0), 0); // Mean == zero
+
+ skew_normal_distribution<RealType> defsd_norm01(0); // Test default (0, sd = 1)
+ BOOST_CHECK_CLOSE(
+ mean(defsd_norm01),
+ static_cast<RealType>(0), 0); // Mean == zero
+
+ skew_normal_distribution<RealType> def_norm01; // Test default (0, sd = 1)
+ BOOST_CHECK_CLOSE(
+ mean(def_norm01),
+ static_cast<RealType>(0), 0); // Mean == zero
+
+ BOOST_CHECK_CLOSE(
+ standard_deviation(def_norm01),
+ static_cast<RealType>(1), 0); //
+
+ BOOST_CHECK_CLOSE(
+ mode(def_norm01),
+ static_cast<RealType>(0), 0); // Mode == zero
+
+
+ // Skew_normal tests with shape != 0.
+ {
+ RealType tol5 = boost::math::tools::epsilon<RealType>() * 5;
+ RealType tol100 = boost::math::tools::epsilon<RealType>() * 100;
+ RealType tol1000 = boost::math::tools::epsilon<RealType>() * 1000;
+
+ skew_normal_distribution<RealType> dist(1.1, 0.02, 0.03);
+
+ BOOST_MATH_STD_USING // ADL of std math lib names.
+
+ // Test values from R = see skew_normal_drv.cpp which included the R code used.
+ {
+ skew_normal_distribution<RealType> dist(1.1, 2.2, -3.3);
+
+ BOOST_CHECK_CLOSE( // mean:
+ mean(dist)
+ , static_cast<RealType>(-0.57990899253985684L), tol1000);
+ BOOST_CHECK_CLOSE( // variance:
+ variance(dist)
+ , static_cast<RealType>(2.017905776783723L), tol100);
+
+ BOOST_CHECK_CLOSE( // skewness:
+ skewness(dist)
+ , static_cast<RealType>(-0.70985454817153848L), tol1000);
+ BOOST_CHECK_CLOSE( // kurtosis:
+ kurtosis(dist)
+ , static_cast<RealType>(3.L + 0.55387526252417818L), tol1000);
+ BOOST_CHECK_CLOSE( // kurtosis excess:
+ kurtosis_excess(dist)
+ , static_cast<RealType>(0.55387526252417818L), tol1000);
+
+ BOOST_CHECK_CLOSE(
+ pdf(dist, static_cast<RealType>(0.4)),
+ static_cast<RealType>(0.29414011015659952L),
+ tol100);
+
+ BOOST_CHECK_CLOSE(
+ cdf(dist, static_cast<RealType>(0.4L)),
+ static_cast<RealType>(0.73391861892787402L),
+ tol1000);
+
+ BOOST_CHECK_CLOSE(
+ quantile(dist, static_cast<RealType>(0.3L)),
+ static_cast<RealType>(-1.180104068086876L),
+ tol1000);
+
+
+ { // mode tests
+
+ skew_normal_distribution<RealType> dist(0, 1, 4.);
+
+ cout << "pdf(dist, 0) = " << pdf(dist, 0) << ", pdf(dist, 0.45) = " << pdf(dist, 0.45) << endl;
+ // BOOST_CHECK_CLOSE(mode(dist), boost::math::constants::root_two<RealType>() / 2, tol5);
+ BOOST_CHECK_CLOSE(mode(dist), 0.695, tol5);
+ }
+
+
+ }
+ {
+ skew_normal_distribution<RealType> dist(1.1, 0.02, 0.03);
+
+ BOOST_CHECK_CLOSE( // mean:
+ mean(dist)
+ , static_cast<RealType>(1.1004785154529559L), tol5);
+ BOOST_CHECK_CLOSE( // variance:
+ variance(dist)
+ , static_cast<RealType>(0.00039977102296128255L), tol100);
+
+ BOOST_CHECK_CLOSE( // skewness:
+ skewness(dist)
+ , static_cast<RealType>(5.8834811259890419e-006L), tol1000);
+ BOOST_CHECK_CLOSE( // kurtosis:
+ kurtosis(dist)
+ , static_cast<RealType>(3.L + 9.2903475812137596e-008L), tol1000);
+ BOOST_CHECK_CLOSE( // kurtosis excess:
+ kurtosis_excess(dist)
+ , static_cast<RealType>(9.2903475812137596e-008L), tol1000 * 10);
+ }
+ {
+ skew_normal_distribution<RealType> dist(10.1, 5., -0.03);
+ BOOST_CHECK_CLOSE( // mean:
+ mean(dist)
+ , static_cast<RealType>(9.9803711367610521L), tol5);
+ BOOST_CHECK_CLOSE( // variance:
+ variance(dist)
+ , static_cast<RealType>(24.985688935080159L), tol100);
+
+ BOOST_CHECK_CLOSE( // skewness:
+ skewness(dist)
+ , static_cast<RealType>(-5.8834811259890411e-006L), tol1000);
+ BOOST_CHECK_CLOSE( // kurtosis:
+ kurtosis(dist)
+ , static_cast<RealType>(3.L + 9.2903475812137596e-008L), tol1000);
+ BOOST_CHECK_CLOSE( // kurtosis excess:
+ kurtosis_excess(dist)
+ , static_cast<RealType>(9.2903475812137596e-008L), tol1000);
+ }
+ {
+ skew_normal_distribution<RealType> dist(-10.1, 5., 30.);
+ BOOST_CHECK_CLOSE( // mean:
+ mean(dist)
+ , static_cast<RealType>(-6.1127916967413842L), tol100);
+ BOOST_CHECK_CLOSE( // variance:
+ variance(dist)
+ , static_cast<RealType>(9.1021699464255477L), tol100);
+
+ BOOST_CHECK_CLOSE( // skewness:
+ skewness(dist)
+ , static_cast<RealType>(0.99072425443686996L), tol1000);
+ BOOST_CHECK_CLOSE( // kurtosis:
+ kurtosis(dist)
+ , static_cast<RealType>(3.L + 0.86388620084060674L), tol1000);
+ BOOST_CHECK_CLOSE( // kurtosis excess:
+ kurtosis_excess(dist)
+ , static_cast<RealType>(0.86388620084060674L), tol1000);
+ }
+
+
+ }
+
+
+} // template <class RealType>void test_spots(RealType)
+
+int test_main(int, char* [])
+{
+
+
+ using boost::math::skew_normal;
+ using boost::math::skew_normal_distribution;
+
+ //int precision = 17; // std::numeric_limits<double::max_digits10;
+ double tolfeweps = numeric_limits<double>::epsilon() * 5;
+ //double tol6decdigits = numeric_limits<float>::epsilon() * 2;
+ // Check that can generate skew_normal distribution using the two convenience methods:
+ boost::math::skew_normal w12(1., 2); // Using typedef.
+ boost::math::skew_normal_distribution<> w01; // Use default unity values for mean and scale.
+ // Note NOT myn01() as the compiler will interpret as a function!
+
+ // Checks on constructors.
+ // Default parameters.
+ BOOST_CHECK_EQUAL(w01.location(), 0);
+ BOOST_CHECK_EQUAL(w01.scale(), 1);
+ BOOST_CHECK_EQUAL(w01.shape(), 0);
+
+ skew_normal_distribution<> w23(2., 3); // Using default RealType double.
+ BOOST_CHECK_EQUAL(w23.scale(), 3);
+ BOOST_CHECK_EQUAL(w23.shape(), 0);
+
+ skew_normal_distribution<> w123(1., 2., 3.); // Using default RealType double.
+ BOOST_CHECK_EQUAL(w123.location(), 1.);
+ BOOST_CHECK_EQUAL(w123.scale(), 2.);
+ BOOST_CHECK_EQUAL(w123.shape(), 3.);
+
+ BOOST_CHECK_CLOSE_FRACTION(mean(w01), static_cast<double>(0), tolfeweps); // Default mean == zero
+ BOOST_CHECK_CLOSE_FRACTION(scale(w01), static_cast<double>(1), tolfeweps); // Default scale == unity
+
+ // Basic sanity-check spot values for all floating-point types..
+ // (Parameter value, arbitrarily zero, only communicates the floating point type).
+ test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 %
+ test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 %
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+ test_spots(0.0L); // Test long double.
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x0582))
+ test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
+#endif
+#else
+ std::cout << "<note>The long double tests have been disabled on this platform "
+ "either because the long double overloads of the usual math functions are "
+ "not available at all, or because they are too inaccurate for these tests "
+ "to pass.</note>" << std::cout;
+#endif
+ /* */
+ return 0;
+} // int test_main(int, char* [])
+
+/*
+
+Output:
+
+
+*/
+
+