276 lines
8.0 KiB
C++
276 lines
8.0 KiB
C++
// 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
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|