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