eval.pass.cpp (23871B)
1 //===----------------------------------------------------------------------===// 2 // 3 // The LLVM Compiler Infrastructure 4 // 5 // This file is dual licensed under the MIT and the University of Illinois Open 6 // Source Licenses. See LICENSE.TXT for details. 7 // 8 //===----------------------------------------------------------------------===// 9 // 10 // REQUIRES: long_tests 11 12 // <random> 13 14 // template<class RealType = double> 15 // class piecewise_constant_distribution 16 17 // template<class _URNG> result_type operator()(_URNG& g); 18 19 #include <random> 20 #include <vector> 21 #include <iterator> 22 #include <numeric> 23 #include <algorithm> // for sort 24 #include <cassert> 25 26 template <class T> 27 inline 28 T 29 sqr(T x) 30 { 31 return x*x; 32 } 33 34 void 35 test1() 36 { 37 typedef std::piecewise_constant_distribution<> D; 38 typedef std::mt19937_64 G; 39 G g; 40 double b[] = {10, 14, 16, 17}; 41 double p[] = {25, 62.5, 12.5}; 42 const size_t Np = sizeof(p) / sizeof(p[0]); 43 D d(b, b+Np+1, p); 44 const int N = 1000000; 45 std::vector<D::result_type> u; 46 for (int i = 0; i < N; ++i) 47 { 48 D::result_type v = d(g); 49 assert(d.min() <= v && v < d.max()); 50 u.push_back(v); 51 } 52 std::vector<double> prob(std::begin(p), std::end(p)); 53 double s = std::accumulate(prob.begin(), prob.end(), 0.0); 54 for (std::size_t i = 0; i < prob.size(); ++i) 55 prob[i] /= s; 56 std::sort(u.begin(), u.end()); 57 for (std::size_t i = 0; i < Np; ++i) 58 { 59 typedef std::vector<D::result_type>::iterator I; 60 I lb = std::lower_bound(u.begin(), u.end(), b[i]); 61 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]); 62 const size_t Ni = ub - lb; 63 if (prob[i] == 0) 64 assert(Ni == 0); 65 else 66 { 67 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01); 68 double mean = std::accumulate(lb, ub, 0.0) / Ni; 69 double var = 0; 70 double skew = 0; 71 double kurtosis = 0; 72 for (I j = lb; j != ub; ++j) 73 { 74 double dbl = (*j - mean); 75 double d2 = sqr(dbl); 76 var += d2; 77 skew += dbl * d2; 78 kurtosis += d2 * d2; 79 } 80 var /= Ni; 81 double dev = std::sqrt(var); 82 skew /= Ni * dev * var; 83 kurtosis /= Ni * var * var; 84 kurtosis -= 3; 85 double x_mean = (b[i+1] + b[i]) / 2; 86 double x_var = sqr(b[i+1] - b[i]) / 12; 87 double x_skew = 0; 88 double x_kurtosis = -6./5; 89 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 90 assert(std::abs((var - x_var) / x_var) < 0.01); 91 assert(std::abs(skew - x_skew) < 0.01); 92 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 93 } 94 } 95 } 96 97 void 98 test2() 99 { 100 typedef std::piecewise_constant_distribution<> D; 101 typedef std::mt19937_64 G; 102 G g; 103 double b[] = {10, 14, 16, 17}; 104 double p[] = {0, 62.5, 12.5}; 105 const size_t Np = sizeof(p) / sizeof(p[0]); 106 D d(b, b+Np+1, p); 107 const int N = 1000000; 108 std::vector<D::result_type> u; 109 for (int i = 0; i < N; ++i) 110 { 111 D::result_type v = d(g); 112 assert(d.min() <= v && v < d.max()); 113 u.push_back(v); 114 } 115 std::vector<double> prob(std::begin(p), std::end(p)); 116 double s = std::accumulate(prob.begin(), prob.end(), 0.0); 117 for (std::size_t i = 0; i < prob.size(); ++i) 118 prob[i] /= s; 119 std::sort(u.begin(), u.end()); 120 for (std::size_t i = 0; i < Np; ++i) 121 { 122 typedef std::vector<D::result_type>::iterator I; 123 I lb = std::lower_bound(u.begin(), u.end(), b[i]); 124 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]); 125 const size_t Ni = ub - lb; 126 if (prob[i] == 0) 127 assert(Ni == 0); 128 else 129 { 130 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01); 131 double mean = std::accumulate(lb, ub, 0.0) / Ni; 132 double var = 0; 133 double skew = 0; 134 double kurtosis = 0; 135 for (I j = lb; j != ub; ++j) 136 { 137 double dbl = (*j - mean); 138 double d2 = sqr(dbl); 139 var += d2; 140 skew += dbl * d2; 141 kurtosis += d2 * d2; 142 } 143 var /= Ni; 144 double dev = std::sqrt(var); 145 skew /= Ni * dev * var; 146 kurtosis /= Ni * var * var; 147 kurtosis -= 3; 148 double x_mean = (b[i+1] + b[i]) / 2; 149 double x_var = sqr(b[i+1] - b[i]) / 12; 150 double x_skew = 0; 151 double x_kurtosis = -6./5; 152 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 153 assert(std::abs((var - x_var) / x_var) < 0.01); 154 assert(std::abs(skew - x_skew) < 0.01); 155 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 156 } 157 } 158 } 159 160 void 161 test3() 162 { 163 typedef std::piecewise_constant_distribution<> D; 164 typedef std::mt19937_64 G; 165 G g; 166 double b[] = {10, 14, 16, 17}; 167 double p[] = {25, 0, 12.5}; 168 const size_t Np = sizeof(p) / sizeof(p[0]); 169 D d(b, b+Np+1, p); 170 const int N = 1000000; 171 std::vector<D::result_type> u; 172 for (int i = 0; i < N; ++i) 173 { 174 D::result_type v = d(g); 175 assert(d.min() <= v && v < d.max()); 176 u.push_back(v); 177 } 178 std::vector<double> prob(std::begin(p), std::end(p)); 179 double s = std::accumulate(prob.begin(), prob.end(), 0.0); 180 for (std::size_t i = 0; i < prob.size(); ++i) 181 prob[i] /= s; 182 std::sort(u.begin(), u.end()); 183 for (std::size_t i = 0; i < Np; ++i) 184 { 185 typedef std::vector<D::result_type>::iterator I; 186 I lb = std::lower_bound(u.begin(), u.end(), b[i]); 187 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]); 188 const size_t Ni = ub - lb; 189 if (prob[i] == 0) 190 assert(Ni == 0); 191 else 192 { 193 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01); 194 double mean = std::accumulate(lb, ub, 0.0) / Ni; 195 double var = 0; 196 double skew = 0; 197 double kurtosis = 0; 198 for (I j = lb; j != ub; ++j) 199 { 200 double dbl = (*j - mean); 201 double d2 = sqr(dbl); 202 var += d2; 203 skew += dbl * d2; 204 kurtosis += d2 * d2; 205 } 206 var /= Ni; 207 double dev = std::sqrt(var); 208 skew /= Ni * dev * var; 209 kurtosis /= Ni * var * var; 210 kurtosis -= 3; 211 double x_mean = (b[i+1] + b[i]) / 2; 212 double x_var = sqr(b[i+1] - b[i]) / 12; 213 double x_skew = 0; 214 double x_kurtosis = -6./5; 215 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 216 assert(std::abs((var - x_var) / x_var) < 0.01); 217 assert(std::abs(skew - x_skew) < 0.01); 218 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 219 } 220 } 221 } 222 223 void 224 test4() 225 { 226 typedef std::piecewise_constant_distribution<> D; 227 typedef std::mt19937_64 G; 228 G g; 229 double b[] = {10, 14, 16, 17}; 230 double p[] = {25, 62.5, 0}; 231 const size_t Np = sizeof(p) / sizeof(p[0]); 232 D d(b, b+Np+1, p); 233 const int N = 1000000; 234 std::vector<D::result_type> u; 235 for (int i = 0; i < N; ++i) 236 { 237 D::result_type v = d(g); 238 assert(d.min() <= v && v < d.max()); 239 u.push_back(v); 240 } 241 std::vector<double> prob(std::begin(p), std::end(p)); 242 double s = std::accumulate(prob.begin(), prob.end(), 0.0); 243 for (std::size_t i = 0; i < prob.size(); ++i) 244 prob[i] /= s; 245 std::sort(u.begin(), u.end()); 246 for (std::size_t i = 0; i < Np; ++i) 247 { 248 typedef std::vector<D::result_type>::iterator I; 249 I lb = std::lower_bound(u.begin(), u.end(), b[i]); 250 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]); 251 const size_t Ni = ub - lb; 252 if (prob[i] == 0) 253 assert(Ni == 0); 254 else 255 { 256 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01); 257 double mean = std::accumulate(lb, ub, 0.0) / Ni; 258 double var = 0; 259 double skew = 0; 260 double kurtosis = 0; 261 for (I j = lb; j != ub; ++j) 262 { 263 double dbl = (*j - mean); 264 double d2 = sqr(dbl); 265 var += d2; 266 skew += dbl * d2; 267 kurtosis += d2 * d2; 268 } 269 var /= Ni; 270 double dev = std::sqrt(var); 271 skew /= Ni * dev * var; 272 kurtosis /= Ni * var * var; 273 kurtosis -= 3; 274 double x_mean = (b[i+1] + b[i]) / 2; 275 double x_var = sqr(b[i+1] - b[i]) / 12; 276 double x_skew = 0; 277 double x_kurtosis = -6./5; 278 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 279 assert(std::abs((var - x_var) / x_var) < 0.01); 280 assert(std::abs(skew - x_skew) < 0.01); 281 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 282 } 283 } 284 } 285 286 void 287 test5() 288 { 289 typedef std::piecewise_constant_distribution<> D; 290 typedef std::mt19937_64 G; 291 G g; 292 double b[] = {10, 14, 16, 17}; 293 double p[] = {25, 0, 0}; 294 const size_t Np = sizeof(p) / sizeof(p[0]); 295 D d(b, b+Np+1, p); 296 const int N = 100000; 297 std::vector<D::result_type> u; 298 for (int i = 0; i < N; ++i) 299 { 300 D::result_type v = d(g); 301 assert(d.min() <= v && v < d.max()); 302 u.push_back(v); 303 } 304 std::vector<double> prob(std::begin(p), std::end(p)); 305 double s = std::accumulate(prob.begin(), prob.end(), 0.0); 306 for (std::size_t i = 0; i < prob.size(); ++i) 307 prob[i] /= s; 308 std::sort(u.begin(), u.end()); 309 for (std::size_t i = 0; i < Np; ++i) 310 { 311 typedef std::vector<D::result_type>::iterator I; 312 I lb = std::lower_bound(u.begin(), u.end(), b[i]); 313 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]); 314 const size_t Ni = ub - lb; 315 if (prob[i] == 0) 316 assert(Ni == 0); 317 else 318 { 319 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01); 320 double mean = std::accumulate(lb, ub, 0.0) / Ni; 321 double var = 0; 322 double skew = 0; 323 double kurtosis = 0; 324 for (I j = lb; j != ub; ++j) 325 { 326 double dbl = (*j - mean); 327 double d2 = sqr(dbl); 328 var += d2; 329 skew += dbl * d2; 330 kurtosis += d2 * d2; 331 } 332 var /= Ni; 333 double dev = std::sqrt(var); 334 skew /= Ni * dev * var; 335 kurtosis /= Ni * var * var; 336 kurtosis -= 3; 337 double x_mean = (b[i+1] + b[i]) / 2; 338 double x_var = sqr(b[i+1] - b[i]) / 12; 339 double x_skew = 0; 340 double x_kurtosis = -6./5; 341 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 342 assert(std::abs((var - x_var) / x_var) < 0.01); 343 assert(std::abs(skew - x_skew) < 0.01); 344 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 345 } 346 } 347 } 348 349 void 350 test6() 351 { 352 typedef std::piecewise_constant_distribution<> D; 353 typedef std::mt19937_64 G; 354 G g; 355 double b[] = {10, 14, 16, 17}; 356 double p[] = {0, 25, 0}; 357 const size_t Np = sizeof(p) / sizeof(p[0]); 358 D d(b, b+Np+1, p); 359 const int N = 100000; 360 std::vector<D::result_type> u; 361 for (int i = 0; i < N; ++i) 362 { 363 D::result_type v = d(g); 364 assert(d.min() <= v && v < d.max()); 365 u.push_back(v); 366 } 367 std::vector<double> prob(std::begin(p), std::end(p)); 368 double s = std::accumulate(prob.begin(), prob.end(), 0.0); 369 for (std::size_t i = 0; i < prob.size(); ++i) 370 prob[i] /= s; 371 std::sort(u.begin(), u.end()); 372 for (std::size_t i = 0; i < Np; ++i) 373 { 374 typedef std::vector<D::result_type>::iterator I; 375 I lb = std::lower_bound(u.begin(), u.end(), b[i]); 376 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]); 377 const size_t Ni = ub - lb; 378 if (prob[i] == 0) 379 assert(Ni == 0); 380 else 381 { 382 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01); 383 double mean = std::accumulate(lb, ub, 0.0) / Ni; 384 double var = 0; 385 double skew = 0; 386 double kurtosis = 0; 387 for (I j = lb; j != ub; ++j) 388 { 389 double dbl = (*j - mean); 390 double d2 = sqr(dbl); 391 var += d2; 392 skew += dbl * d2; 393 kurtosis += d2 * d2; 394 } 395 var /= Ni; 396 double dev = std::sqrt(var); 397 skew /= Ni * dev * var; 398 kurtosis /= Ni * var * var; 399 kurtosis -= 3; 400 double x_mean = (b[i+1] + b[i]) / 2; 401 double x_var = sqr(b[i+1] - b[i]) / 12; 402 double x_skew = 0; 403 double x_kurtosis = -6./5; 404 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 405 assert(std::abs((var - x_var) / x_var) < 0.01); 406 assert(std::abs(skew - x_skew) < 0.01); 407 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 408 } 409 } 410 } 411 412 void 413 test7() 414 { 415 typedef std::piecewise_constant_distribution<> D; 416 typedef std::mt19937_64 G; 417 G g; 418 double b[] = {10, 14, 16, 17}; 419 double p[] = {0, 0, 1}; 420 const size_t Np = sizeof(p) / sizeof(p[0]); 421 D d(b, b+Np+1, p); 422 const int N = 100000; 423 std::vector<D::result_type> u; 424 for (int i = 0; i < N; ++i) 425 { 426 D::result_type v = d(g); 427 assert(d.min() <= v && v < d.max()); 428 u.push_back(v); 429 } 430 std::vector<double> prob(std::begin(p), std::end(p)); 431 double s = std::accumulate(prob.begin(), prob.end(), 0.0); 432 for (std::size_t i = 0; i < prob.size(); ++i) 433 prob[i] /= s; 434 std::sort(u.begin(), u.end()); 435 for (std::size_t i = 0; i < Np; ++i) 436 { 437 typedef std::vector<D::result_type>::iterator I; 438 I lb = std::lower_bound(u.begin(), u.end(), b[i]); 439 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]); 440 const size_t Ni = ub - lb; 441 if (prob[i] == 0) 442 assert(Ni == 0); 443 else 444 { 445 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01); 446 double mean = std::accumulate(lb, ub, 0.0) / Ni; 447 double var = 0; 448 double skew = 0; 449 double kurtosis = 0; 450 for (I j = lb; j != ub; ++j) 451 { 452 double dbl = (*j - mean); 453 double d2 = sqr(dbl); 454 var += d2; 455 skew += dbl * d2; 456 kurtosis += d2 * d2; 457 } 458 var /= Ni; 459 double dev = std::sqrt(var); 460 skew /= Ni * dev * var; 461 kurtosis /= Ni * var * var; 462 kurtosis -= 3; 463 double x_mean = (b[i+1] + b[i]) / 2; 464 double x_var = sqr(b[i+1] - b[i]) / 12; 465 double x_skew = 0; 466 double x_kurtosis = -6./5; 467 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 468 assert(std::abs((var - x_var) / x_var) < 0.01); 469 assert(std::abs(skew - x_skew) < 0.01); 470 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 471 } 472 } 473 } 474 475 void 476 test8() 477 { 478 typedef std::piecewise_constant_distribution<> D; 479 typedef std::mt19937_64 G; 480 G g; 481 double b[] = {10, 14, 16}; 482 double p[] = {75, 25}; 483 const size_t Np = sizeof(p) / sizeof(p[0]); 484 D d(b, b+Np+1, p); 485 const int N = 100000; 486 std::vector<D::result_type> u; 487 for (int i = 0; i < N; ++i) 488 { 489 D::result_type v = d(g); 490 assert(d.min() <= v && v < d.max()); 491 u.push_back(v); 492 } 493 std::vector<double> prob(std::begin(p), std::end(p)); 494 double s = std::accumulate(prob.begin(), prob.end(), 0.0); 495 for (std::size_t i = 0; i < prob.size(); ++i) 496 prob[i] /= s; 497 std::sort(u.begin(), u.end()); 498 for (std::size_t i = 0; i < Np; ++i) 499 { 500 typedef std::vector<D::result_type>::iterator I; 501 I lb = std::lower_bound(u.begin(), u.end(), b[i]); 502 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]); 503 const size_t Ni = ub - lb; 504 if (prob[i] == 0) 505 assert(Ni == 0); 506 else 507 { 508 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01); 509 double mean = std::accumulate(lb, ub, 0.0) / Ni; 510 double var = 0; 511 double skew = 0; 512 double kurtosis = 0; 513 for (I j = lb; j != ub; ++j) 514 { 515 double dbl = (*j - mean); 516 double d2 = sqr(dbl); 517 var += d2; 518 skew += dbl * d2; 519 kurtosis += d2 * d2; 520 } 521 var /= Ni; 522 double dev = std::sqrt(var); 523 skew /= Ni * dev * var; 524 kurtosis /= Ni * var * var; 525 kurtosis -= 3; 526 double x_mean = (b[i+1] + b[i]) / 2; 527 double x_var = sqr(b[i+1] - b[i]) / 12; 528 double x_skew = 0; 529 double x_kurtosis = -6./5; 530 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 531 assert(std::abs((var - x_var) / x_var) < 0.01); 532 assert(std::abs(skew - x_skew) < 0.01); 533 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 534 } 535 } 536 } 537 538 void 539 test9() 540 { 541 typedef std::piecewise_constant_distribution<> D; 542 typedef std::mt19937_64 G; 543 G g; 544 double b[] = {10, 14, 16}; 545 double p[] = {0, 25}; 546 const size_t Np = sizeof(p) / sizeof(p[0]); 547 D d(b, b+Np+1, p); 548 const int N = 100000; 549 std::vector<D::result_type> u; 550 for (int i = 0; i < N; ++i) 551 { 552 D::result_type v = d(g); 553 assert(d.min() <= v && v < d.max()); 554 u.push_back(v); 555 } 556 std::vector<double> prob(std::begin(p), std::end(p)); 557 double s = std::accumulate(prob.begin(), prob.end(), 0.0); 558 for (std::size_t i = 0; i < prob.size(); ++i) 559 prob[i] /= s; 560 std::sort(u.begin(), u.end()); 561 for (std::size_t i = 0; i < Np; ++i) 562 { 563 typedef std::vector<D::result_type>::iterator I; 564 I lb = std::lower_bound(u.begin(), u.end(), b[i]); 565 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]); 566 const size_t Ni = ub - lb; 567 if (prob[i] == 0) 568 assert(Ni == 0); 569 else 570 { 571 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01); 572 double mean = std::accumulate(lb, ub, 0.0) / Ni; 573 double var = 0; 574 double skew = 0; 575 double kurtosis = 0; 576 for (I j = lb; j != ub; ++j) 577 { 578 double dbl = (*j - mean); 579 double d2 = sqr(dbl); 580 var += d2; 581 skew += dbl * d2; 582 kurtosis += d2 * d2; 583 } 584 var /= Ni; 585 double dev = std::sqrt(var); 586 skew /= Ni * dev * var; 587 kurtosis /= Ni * var * var; 588 kurtosis -= 3; 589 double x_mean = (b[i+1] + b[i]) / 2; 590 double x_var = sqr(b[i+1] - b[i]) / 12; 591 double x_skew = 0; 592 double x_kurtosis = -6./5; 593 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 594 assert(std::abs((var - x_var) / x_var) < 0.01); 595 assert(std::abs(skew - x_skew) < 0.01); 596 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 597 } 598 } 599 } 600 601 void 602 test10() 603 { 604 typedef std::piecewise_constant_distribution<> D; 605 typedef std::mt19937_64 G; 606 G g; 607 double b[] = {10, 14, 16}; 608 double p[] = {1, 0}; 609 const size_t Np = sizeof(p) / sizeof(p[0]); 610 D d(b, b+Np+1, p); 611 const int N = 100000; 612 std::vector<D::result_type> u; 613 for (int i = 0; i < N; ++i) 614 { 615 D::result_type v = d(g); 616 assert(d.min() <= v && v < d.max()); 617 u.push_back(v); 618 } 619 std::vector<double> prob(std::begin(p), std::end(p)); 620 double s = std::accumulate(prob.begin(), prob.end(), 0.0); 621 for (std::size_t i = 0; i < prob.size(); ++i) 622 prob[i] /= s; 623 std::sort(u.begin(), u.end()); 624 for (std::size_t i = 0; i < Np; ++i) 625 { 626 typedef std::vector<D::result_type>::iterator I; 627 I lb = std::lower_bound(u.begin(), u.end(), b[i]); 628 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]); 629 const size_t Ni = ub - lb; 630 if (prob[i] == 0) 631 assert(Ni == 0); 632 else 633 { 634 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01); 635 double mean = std::accumulate(lb, ub, 0.0) / Ni; 636 double var = 0; 637 double skew = 0; 638 double kurtosis = 0; 639 for (I j = lb; j != ub; ++j) 640 { 641 double dbl = (*j - mean); 642 double d2 = sqr(dbl); 643 var += d2; 644 skew += dbl * d2; 645 kurtosis += d2 * d2; 646 } 647 var /= Ni; 648 double dev = std::sqrt(var); 649 skew /= Ni * dev * var; 650 kurtosis /= Ni * var * var; 651 kurtosis -= 3; 652 double x_mean = (b[i+1] + b[i]) / 2; 653 double x_var = sqr(b[i+1] - b[i]) / 12; 654 double x_skew = 0; 655 double x_kurtosis = -6./5; 656 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 657 assert(std::abs((var - x_var) / x_var) < 0.01); 658 assert(std::abs(skew - x_skew) < 0.01); 659 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 660 } 661 } 662 } 663 664 void 665 test11() 666 { 667 typedef std::piecewise_constant_distribution<> D; 668 typedef std::mt19937_64 G; 669 G g; 670 double b[] = {10, 14}; 671 double p[] = {1}; 672 const size_t Np = sizeof(p) / sizeof(p[0]); 673 D d(b, b+Np+1, p); 674 const int N = 100000; 675 std::vector<D::result_type> u; 676 for (int i = 0; i < N; ++i) 677 { 678 D::result_type v = d(g); 679 assert(d.min() <= v && v < d.max()); 680 u.push_back(v); 681 } 682 std::vector<double> prob(std::begin(p), std::end(p)); 683 double s = std::accumulate(prob.begin(), prob.end(), 0.0); 684 for (std::size_t i = 0; i < prob.size(); ++i) 685 prob[i] /= s; 686 std::sort(u.begin(), u.end()); 687 for (std::size_t i = 0; i < Np; ++i) 688 { 689 typedef std::vector<D::result_type>::iterator I; 690 I lb = std::lower_bound(u.begin(), u.end(), b[i]); 691 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]); 692 const size_t Ni = ub - lb; 693 if (prob[i] == 0) 694 assert(Ni == 0); 695 else 696 { 697 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01); 698 double mean = std::accumulate(lb, ub, 0.0) / Ni; 699 double var = 0; 700 double skew = 0; 701 double kurtosis = 0; 702 for (I j = lb; j != ub; ++j) 703 { 704 double dbl = (*j - mean); 705 double d2 = sqr(dbl); 706 var += d2; 707 skew += dbl * d2; 708 kurtosis += d2 * d2; 709 } 710 var /= Ni; 711 double dev = std::sqrt(var); 712 skew /= Ni * dev * var; 713 kurtosis /= Ni * var * var; 714 kurtosis -= 3; 715 double x_mean = (b[i+1] + b[i]) / 2; 716 double x_var = sqr(b[i+1] - b[i]) / 12; 717 double x_skew = 0; 718 double x_kurtosis = -6./5; 719 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 720 assert(std::abs((var - x_var) / x_var) < 0.01); 721 assert(std::abs(skew - x_skew) < 0.01); 722 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 723 } 724 } 725 } 726 727 int main() 728 { 729 test1(); 730 test2(); 731 test3(); 732 test4(); 733 test5(); 734 test6(); 735 test7(); 736 test8(); 737 test9(); 738 test10(); 739 test11(); 740 }