projects / ceph.git / blame
| Commit | Line | Data |
|---|---|---|
| 7c673cae FG |
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 | ||
| 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 | ||
| 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 |