]>
git.proxmox.com Git - ceph.git/blob - ceph/src/boost/libs/math/example/skew_normal_example.cpp
1 // Copyright Benjamin Sobotta 2012
3 // Use, modification and distribution are subject to the
4 // Boost Software License, Version 1.0. (See accompanying file
5 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
9 # pragma warning (disable : 4512) // assignment operator could not be generated
10 # pragma warning (disable : 4127) // conditional expression is constant.
13 #include <boost/math/distributions/skew_normal.hpp>
14 using boost::math::skew_normal_distribution
;
15 using boost::math::skew_normal
;
20 void check(const double loc
, const double sc
, const double sh
,
21 const double * const cumulants
, const std::pair
<double, double> qu
,
22 const double x
, const double tpdf
, const double tcdf
)
24 using namespace boost::math
;
26 skew_normal
D(loc
, sc
, sh
);
28 const double sk
= cumulants
[2] / (std::pow(cumulants
[1], 1.5));
29 const double kt
= cumulants
[3] / (cumulants
[1] * cumulants
[1]);
31 // checks against tabulated values
32 std::cout
<< "mean: table=" << cumulants
[0] << "\tcompute=" << mean(D
) << "\tdiff=" << fabs(cumulants
[0]-mean(D
)) << std::endl
;
33 std::cout
<< "var: table=" << cumulants
[1] << "\tcompute=" << variance(D
) << "\tdiff=" << fabs(cumulants
[1]-variance(D
)) << std::endl
;
34 std::cout
<< "skew: table=" << sk
<< "\tcompute=" << skewness(D
) << "\tdiff=" << fabs(sk
-skewness(D
)) << std::endl
;
35 std::cout
<< "kur.: table=" << kt
<< "\tcompute=" << kurtosis_excess(D
) << "\tdiff=" << fabs(kt
-kurtosis_excess(D
)) << std::endl
;
36 std::cout
<< "mode: table=" << "N/A" << "\tcompute=" << mode(D
) << "\tdiff=" << "N/A" << std::endl
;
38 const double q
= quantile(D
, qu
.first
);
39 const double cq
= quantile(complement(D
, qu
.first
));
41 std::cout
<< "quantile(" << qu
.first
<< "): table=" << qu
.second
<< "\tcompute=" << q
<< "\tdiff=" << fabs(qu
.second
-q
) << std::endl
;
44 std::cout
<< "cdf(quantile)=" << cdf(D
, q
) << "\tp=" << qu
.first
<< "\tdiff=" << fabs(qu
.first
-cdf(D
, q
)) << std::endl
;
45 std::cout
<< "ccdf(cquantile)=" << cdf(complement(D
,cq
)) << "\tp=" << qu
.first
<< "\tdiff=" << fabs(qu
.first
-cdf(complement(D
,cq
))) << std::endl
;
48 std::cout
<< "pdf(" << x
<< "): table=" << tpdf
<< "\tcompute=" << pdf(D
,x
) << "\tdiff=" << fabs(tpdf
-pdf(D
,x
)) << std::endl
;
49 std::cout
<< "cdf(" << x
<< "): table=" << tcdf
<< "\tcompute=" << cdf(D
,x
) << "\tdiff=" << fabs(tcdf
-cdf(D
,x
)) << std::endl
;
50 std::cout
<< "================================\n";
55 using namespace boost::math
;
57 double sc
= 0.0,loc
,sh
,x
,dsn
,qsn
,psn
,p
;
58 std::cout
<< std::setprecision(20);
64 > install.packages("sn")
65 Warning in install.packages("sn") :
66 'lib = "/usr/lib64/R/library"' is not writable
67 Would you like to create a personal library
68 '~/R/x86_64-unknown-linux-gnu-library/2.12'
69 to install packages into? (y/n) y
70 --- Please select a CRAN mirror for use in this session ---
71 Loading Tcl/Tk interface ... done
72 also installing the dependency mnormt
74 trying URL 'http://mirrors.dotsrc.org/cran/src/contrib/mnormt_1.4-5.tar.gz'
75 Content type 'application/x-gzip' length 34049 bytes (33 Kb)
77 ==================================================
80 trying URL 'http://mirrors.dotsrc.org/cran/src/contrib/sn_0.4-17.tar.gz'
81 Content type 'application/x-gzip' length 65451 bytes (63 Kb)
83 ==================================================
91 > sn.cumulants(1.1, 2.2, -3.3)
92 [1] -0.5799089925398568 2.0179057767837230 -2.0347951542374196
93 [4] 2.2553488991015072
94 > qsn(0.3, 1.1, 2.2, -3.3)
95 [1] -1.180104068086876
96 > psn(0.4, 1.1, 2.2, -3.3)
98 > dsn(0.4, 1.1, 2.2, -3.3)
99 [1] 0.2941401101565995
104 loc
= 1.1; sc
= 2.2; sh
= -3.3;
105 std::cout
<< "location: " << loc
<< "\tscale: " << sc
<< "\tshape: " << sh
<< std::endl
;
106 cumulants
[0] = -0.5799089925398568;
107 cumulants
[1] = 2.0179057767837230;
108 cumulants
[2] = -2.0347951542374196;
109 cumulants
[3] = 2.2553488991015072;
112 qsn
= -1.180104068086876;
113 psn
= 0.733918618927874;
114 dsn
= 0.2941401101565995;
116 check(loc
, sc
, sh
, cumulants
, std::make_pair(p
,qsn
), x
, dsn
, psn
);
120 > sn.cumulants(1.1, .02, .03)
121 [1] 1.1004785154529559e+00 3.9977102296128255e-04 4.7027439329779991e-11
122 [4] 1.4847542790693825e-14
123 > qsn(0.01, 1.1, .02, .03)
124 [1] 1.053964962950150
125 > psn(1.3, 1.1, .02, .03)
127 > dsn(1.3, 1.1, .02, .03)
128 [1] 4.754580380601393e-21
133 loc
= 1.1; sc
= .02; sh
= .03;
134 std::cout
<< "location: " << loc
<< "\tscale: " << sc
<< "\tshape: " << sh
<< std::endl
;
135 cumulants
[0] = 1.1004785154529559e+00;
136 cumulants
[1] = 3.9977102296128255e-04;
137 cumulants
[2] = 4.7027439329779991e-11;
138 cumulants
[3] = 1.4847542790693825e-14;
141 qsn
= 1.053964962950150;
143 dsn
= 4.754580380601393e-21;
145 check(loc
, sc
, sh
, cumulants
, std::make_pair(p
,qsn
), x
, dsn
, psn
);
149 > sn.cumulants(10.1, 5, -.03)
150 [1] 9.980371136761052e+00 2.498568893508016e+01 -7.348037395278123e-04
151 [4] 5.799821402614775e-05
152 > qsn(.8, 10.1, 5, -.03)
153 [1] 14.18727407491953
154 > psn(-1.3, 10.1, 5, -.03)
155 [1] 0.01201290665838824
156 > dsn(-1.3, 10.1, 5, -.03)
157 [1] 0.006254346348897927
163 loc
= 10.1; sc
= 5; sh
= -.03;
164 std::cout
<< "location: " << loc
<< "\tscale: " << sc
<< "\tshape: " << sh
<< std::endl
;
165 cumulants
[0] = 9.980371136761052e+00;
166 cumulants
[1] = 2.498568893508016e+01;
167 cumulants
[2] = -7.348037395278123e-04;
168 cumulants
[3] = 5.799821402614775e-05;
171 qsn
= 14.18727407491953;
172 psn
= 0.01201290665838824;
173 dsn
= 0.006254346348897927;
175 check(loc
, sc
, sh
, cumulants
, std::make_pair(p
,qsn
), x
, dsn
, psn
);
180 > sn.cumulants(-10.1, 5, 30)
181 [1] -6.112791696741384 9.102169946425548 27.206345266148194 71.572537838642916
182 > qsn(.8, -10.1, 5, 30)
184 > psn(-1.3, -10.1, 5, 30)
185 [1] 0.921592193425035
186 > dsn(-1.3, -10.1, 5, 30)
187 [1] 0.0339105445232089
192 loc
= -10.1; sc
= 5; sh
= 30;
193 std::cout
<< "location: " << loc
<< "\tscale: " << sc
<< "\tshape: " << sh
<< std::endl
;
194 cumulants
[0] = -6.112791696741384;
195 cumulants
[1] = 9.102169946425548;
196 cumulants
[2] = 27.206345266148194;
197 cumulants
[3] = 71.572537838642916;
200 qsn
= -3.692242172277;
201 psn
= 0.921592193425035;
202 dsn
= 0.0339105445232089;
204 check(loc
, sc
, sh
, cumulants
, std::make_pair(p
,qsn
), x
, dsn
, psn
);
209 > sn.cumulants(0,1,5)
210 [1] 0.7823901817554269 0.3878656034927102 0.2055576317962637 0.1061119471655128
212 [1] 0.674471117502844
214 [1] 0.0002731513884140924
216 [1] 0.00437241570403263
221 loc
= 0; sc
= 1; sh
= 5;
222 std::cout
<< "location: " << loc
<< "\tscale: " << sc
<< "\tshape: " << sh
<< std::endl
;
223 cumulants
[0] = 0.7823901817554269;
224 cumulants
[1] = 0.3878656034927102;
225 cumulants
[2] = 0.2055576317962637;
226 cumulants
[3] = 0.1061119471655128;
229 qsn
= 0.674471117502844;
230 psn
= 0.0002731513884140924;
231 dsn
= 0.00437241570403263;
233 check(loc
, sc
, sh
, cumulants
, std::make_pair(p
,qsn
), x
, dsn
, psn
);
237 > sn.cumulants(0,1,1e5)
238 [1] 0.7978845607629713 0.3633802276960805 0.2180136141122883 0.1147706820312645
240 [1] 0.6744897501960818
249 loc
= 0; sc
= 1; sh
= 1e5
;
250 std::cout
<< "location: " << loc
<< "\tscale: " << sc
<< "\tshape: " << sh
<< std::endl
;
251 cumulants
[0] = 0.7978845607629713;
252 cumulants
[1] = 0.3633802276960805;
253 cumulants
[2] = 0.2180136141122883;
254 cumulants
[3] = 0.1147706820312645;
257 qsn
= 0.6744897501960818;
261 check(loc
, sc
, sh
, cumulants
, std::make_pair(p
,qsn
), x
, dsn
, psn
);