// 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 : 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[1]-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;
std::cout << "mode: table=" << "N/A" << "\tcompute=" << mode(D) << "\tdiff=" << "N/A" << 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;
std::cout << "location: " << loc << "\tscale: " << sc << "\tshape: " << sh << std::endl;
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;
std::cout << "location: " << loc << "\tscale: " << sc << "\tshape: " << sh << std::endl;
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;
std::cout << "location: " << loc << "\tscale: " << sc << "\tshape: " << sh << std::endl;
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;
std::cout << "location: " << loc << "\tscale: " << sc << "\tshape: " << sh << std::endl;
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);
/* R:
> sn.cumulants(0,1,5)
[1] 0.7823901817554269 0.3878656034927102 0.2055576317962637 0.1061119471655128
> qsn(0.5,0,1,5)
[1] 0.674471117502844
> psn(-0.5, 0,1,5)
[1] 0.0002731513884140924
> dsn(-0.5, 0,1,5)
[1] 0.00437241570403263
*/
// 5th
loc = 0; sc = 1; sh = 5;
std::cout << "location: " << loc << "\tscale: " << sc << "\tshape: " << sh << std::endl;
cumulants[0] = 0.7823901817554269;
cumulants[1] = 0.3878656034927102;
cumulants[2] = 0.2055576317962637;
cumulants[3] = 0.1061119471655128;
x = -0.5;
p = 0.5;
qsn = 0.674471117502844;
psn = 0.0002731513884140924;
dsn = 0.00437241570403263;
check(loc, sc, sh, cumulants, std::make_pair(p,qsn), x, dsn, psn);
/* R:
> sn.cumulants(0,1,1e5)
[1] 0.7978845607629713 0.3633802276960805 0.2180136141122883 0.1147706820312645
> qsn(0.5,0,1,1e5)
[1] 0.6744897501960818
> psn(-0.5, 0,1,1e5)
[1] 0
> dsn(-0.5, 0,1,1e5)
[1] 0
*/
// 6th
loc = 0; sc = 1; sh = 1e5;
std::cout << "location: " << loc << "\tscale: " << sc << "\tshape: " << sh << std::endl;
cumulants[0] = 0.7978845607629713;
cumulants[1] = 0.3633802276960805;
cumulants[2] = 0.2180136141122883;
cumulants[3] = 0.1147706820312645;
x = -0.5;
p = 0.5;
qsn = 0.6744897501960818;
psn = 0.;
dsn = 0.;
check(loc, sc, sh, cumulants, std::make_pair(p,qsn), x, dsn, psn);
return 0;
}
// EOF
