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