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