1 // Copyright Benjamin Sobotta 2012
2
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)
6
7
8 #ifdef _MSC_VER
9 # pragma warning (disable : 4512) // assignment operator could not be generated
10 # pragma warning (disable : 4127) // conditional expression is constant.
11 #endif
12
13 #include <boost/math/distributions/skew_normal.hpp>
14 using boost::math::skew_normal_distribution;
15 using boost::math::skew_normal;
16 #include <iostream>
17 #include <cmath>
18 #include <utility>
19
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)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)
23 {
24 using namespace boost::math;
25
26 skew_normal D(loc, sc, sh);
27
28 const double sk = cumulants[2] / (std::pow(cumulants[1], 1.5));
29 const double kt = cumulants[3] / (cumulants[1] * cumulants[1]);
30
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;
37
38 const double q = quantile(D, qu.first);
39 const double cq = quantile(complement(D, qu.first));
40
41 std::cout << "quantile(" << qu.first << "): table=" << qu.second << "\tcompute=" << q << "\tdiff=" << fabs(qu.second-q) << std::endl;
42
43 // consistency
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;
46
47 // PDF & CDF
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";
51 }
52
main()53 int main()
54 {
55 using namespace boost::math;
56
57 double sc = 0.0,loc,sh,x,dsn,qsn,psn,p;
58 std::cout << std::setprecision(20);
59
60 double cumulants[4];
61
62
63 /* R:
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
73
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)
76 opened URL
77 ==================================================
78 downloaded 33 Kb
79
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)
82 opened URL
83 ==================================================
84 downloaded 63 Kb
85
86
87 > library(sn)
88 > options(digits=22)
89
90
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)
97 [1] 0.733918618927874
98 > dsn(0.4, 1.1, 2.2, -3.3)
99 [1] 0.2941401101565995
100
101 */
102
103 //1 st
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;
110 x = 0.4;
111 p = 0.3;
112 qsn = -1.180104068086876;
113 psn = 0.733918618927874;
114 dsn = 0.2941401101565995;
115
116 check(loc, sc, sh, cumulants, std::make_pair(p,qsn), x, dsn, psn);
117
118 /* R:
119
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)
126 [1] 1
127 > dsn(1.3, 1.1, .02, .03)
128 [1] 4.754580380601393e-21
129
130 */
131
132 // 2nd
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;
139 x = 1.3;
140 p = 0.01;
141 qsn = 1.053964962950150;
142 psn = 1;
143 dsn = 4.754580380601393e-21;
144
145 check(loc, sc, sh, cumulants, std::make_pair(p,qsn), x, dsn, psn);
146
147 /* R:
148
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
158
159
160 */
161
162 // 3rd
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;
169 x = -1.3;
170 p = 0.8;
171 qsn = 14.18727407491953;
172 psn = 0.01201290665838824;
173 dsn = 0.006254346348897927;
174
175 check(loc, sc, sh, cumulants, std::make_pair(p,qsn), x, dsn, psn);
176
177
178 /* R:
179
180 > sn.cumulants(-10.1, 5, 30)
181 [1] -6.112791696741384 9.102169946425548 27.206345266148194 71.572537838642916
182 > qsn(.8, -10.1, 5, 30)
183 [1] -3.692242172277
184 > psn(-1.3, -10.1, 5, 30)
185 [1] 0.921592193425035
186 > dsn(-1.3, -10.1, 5, 30)
187 [1] 0.0339105445232089
188
189 */
190
191 // 4th
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;
198 x = -1.3;
199 p = 0.8;
200 qsn = -3.692242172277;
201 psn = 0.921592193425035;
202 dsn = 0.0339105445232089;
203
204 check(loc, sc, sh, cumulants, std::make_pair(p,qsn), x, dsn, psn);
205
206
207 /* R:
208
209 > sn.cumulants(0,1,5)
210 [1] 0.7823901817554269 0.3878656034927102 0.2055576317962637 0.1061119471655128
211 > qsn(0.5,0,1,5)
212 [1] 0.674471117502844
213 > psn(-0.5, 0,1,5)
214 [1] 0.0002731513884140924
215 > dsn(-0.5, 0,1,5)
216 [1] 0.00437241570403263
217
218 */
219
220 // 5th
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;
227 x = -0.5;
228 p = 0.5;
229 qsn = 0.674471117502844;
230 psn = 0.0002731513884140924;
231 dsn = 0.00437241570403263;
232
233 check(loc, sc, sh, cumulants, std::make_pair(p,qsn), x, dsn, psn);
234
235 /* R:
236
237 > sn.cumulants(0,1,1e5)
238 [1] 0.7978845607629713 0.3633802276960805 0.2180136141122883 0.1147706820312645
239 > qsn(0.5,0,1,1e5)
240 [1] 0.6744897501960818
241 > psn(-0.5, 0,1,1e5)
242 [1] 0
243 > dsn(-0.5, 0,1,1e5)
244 [1] 0
245
246 */
247
248 // 6th
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;
255 x = -0.5;
256 p = 0.5;
257 qsn = 0.6744897501960818;
258 psn = 0.;
259 dsn = 0.;
260
261 check(loc, sc, sh, cumulants, std::make_pair(p,qsn), x, dsn, psn);
262
263 return 0;
264 }
265
266
267 // EOF
268
269
270
271
272
273
274
275
276