| // 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 |
| |
| |
| |
| |
| |
| |
| |
| |
