745 lines · cpp
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-exception6//7//===----------------------------------------------------------------------===//8//9// REQUIRES: long_tests10 11// <random>12 13// template<class RealType = double>14// class piecewise_constant_distribution15 16// template<class _URNG> result_type operator()(_URNG& g);17 18#include <random>19#include <algorithm> // for sort20#include <cassert>21#include <cmath>22#include <iterator>23#include <numeric>24#include <vector>25 26#include "test_macros.h"27 28template <class T>29inline30T31sqr(T x)32{33 return x*x;34}35 36void37test1()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 else68 {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 99void100test2()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 else131 {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 162void163test3()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 else194 {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 225void226test4()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 else257 {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 288void289test5()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 else320 {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 351void352test6()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 else383 {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 414void415test7()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 else446 {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 477void478test8()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 else509 {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 540void541test9()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 else572 {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 603void604test10()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 else635 {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 666void667test11()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 else698 {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 729int 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