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