test/std/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval.pass.cpp

   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 IntType = int>
  15 // class binomial_distribution
  16
  17 // template<class _URNG> result_type operator()(_URNG& g);
  18
  19 #include <random>
  20 #include <numeric>
  21 #include <vector>
  22 #include <cassert>
  23
  24 template <class T>
  25 inline
  26 T
  27 sqr(T x)
  28 {
  29     return x * x;
  30 }
  31
  32 void
  33 test1()
  34 {
  35     typedef std::binomial_distribution<> D;
  36     typedef std::mt19937_64 G;
  37     G g;
  38     D d(5, .75);
  39     const int N = 1000000;
  40     std::vector<D::result_type> u;
  41     for (int i = 0; i < N; ++i)
  42     {
  43         D::result_type v = d(g);
  44         assert(d.min() <= v && v <= d.max());
  45         u.push_back(v);
  46     }
  47     double mean = std::accumulate(u.begin(), u.end(),
  48                                           double(0)) / u.size();
  49     double var = 0;
  50     double skew = 0;
  51     double kurtosis = 0;
  52     for (unsigned i = 0; i < u.size(); ++i)
  53     {
  54         double dbl = (u[i] - mean);
  55         double d2 = sqr(dbl);
  56         var += d2;
  57         skew += dbl * d2;
  58         kurtosis += d2 * d2;
  59     }
  60     var /= u.size();
  61     double dev = std::sqrt(var);
  62     skew /= u.size() * dev * var;
  63     kurtosis /= u.size() * var * var;
  64     kurtosis -= 3;
  65     double x_mean = d.t() * d.p();
  66     double x_var = x_mean*(1-d.p());
  67     double x_skew = (1-2*d.p()) / std::sqrt(x_var);
  68     double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
  69     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
  70     assert(std::abs((var - x_var) / x_var) < 0.01);
  71     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
  72     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.04);
  73 }
  74
  75 void
  76 test2()
  77 {
  78     typedef std::binomial_distribution<> D;
  79     typedef std::mt19937 G;
  80     G g;
  81     D d(30, .03125);
  82     const int N = 100000;
  83     std::vector<D::result_type> u;
  84     for (int i = 0; i < N; ++i)
  85     {
  86         D::result_type v = d(g);
  87         assert(d.min() <= v && v <= d.max());
  88         u.push_back(v);
  89     }
  90     double mean = std::accumulate(u.begin(), u.end(),
  91                                           double(0)) / u.size();
  92     double var = 0;
  93     double skew = 0;
  94     double kurtosis = 0;
  95     for (unsigned i = 0; i < u.size(); ++i)
  96     {
  97         double dbl = (u[i] - mean);
  98         double d2 = sqr(dbl);
  99         var += d2;
 100         skew += dbl * d2;
 101         kurtosis += d2 * d2;
 102     }
 103     var /= u.size();
 104     double dev = std::sqrt(var);
 105     skew /= u.size() * dev * var;
 106     kurtosis /= u.size() * var * var;
 107     kurtosis -= 3;
 108     double x_mean = d.t() * d.p();
 109     double x_var = x_mean*(1-d.p());
 110     double x_skew = (1-2*d.p()) / std::sqrt(x_var);
 111     double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
 112     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 113     assert(std::abs((var - x_var) / x_var) < 0.01);
 114     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
 115     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
 116 }
 117
 118 void
 119 test3()
 120 {
 121     typedef std::binomial_distribution<> D;
 122     typedef std::mt19937 G;
 123     G g;
 124     D d(40, .25);
 125     const int N = 100000;
 126     std::vector<D::result_type> u;
 127     for (int i = 0; i < N; ++i)
 128     {
 129         D::result_type v = d(g);
 130         assert(d.min() <= v && v <= d.max());
 131         u.push_back(v);
 132     }
 133     double mean = std::accumulate(u.begin(), u.end(),
 134                                           double(0)) / u.size();
 135     double var = 0;
 136     double skew = 0;
 137     double kurtosis = 0;
 138     for (unsigned i = 0; i < u.size(); ++i)
 139     {
 140         double dbl = (u[i] - mean);
 141         double d2 = sqr(dbl);
 142         var += d2;
 143         skew += dbl * d2;
 144         kurtosis += d2 * d2;
 145     }
 146     var /= u.size();
 147     double dev = std::sqrt(var);
 148     skew /= u.size() * dev * var;
 149     kurtosis /= u.size() * var * var;
 150     kurtosis -= 3;
 151     double x_mean = d.t() * d.p();
 152     double x_var = x_mean*(1-d.p());
 153     double x_skew = (1-2*d.p()) / std::sqrt(x_var);
 154     double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
 155     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 156     assert(std::abs((var - x_var) / x_var) < 0.01);
 157     assert(std::abs((skew - x_skew) / x_skew) < 0.03);
 158     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.3);
 159 }
 160
 161 void
 162 test4()
 163 {
 164     typedef std::binomial_distribution<> D;
 165     typedef std::mt19937 G;
 166     G g;
 167     D d(40, 0);
 168     const int N = 100000;
 169     std::vector<D::result_type> u;
 170     for (int i = 0; i < N; ++i)
 171     {
 172         D::result_type v = d(g);
 173         assert(d.min() <= v && v <= d.max());
 174         u.push_back(v);
 175     }
 176     double mean = std::accumulate(u.begin(), u.end(),
 177                                           double(0)) / u.size();
 178     double var = 0;
 179     double skew = 0;
 180     double kurtosis = 0;
 181     for (unsigned i = 0; i < u.size(); ++i)
 182     {
 183         double dbl = (u[i] - mean);
 184         double d2 = sqr(dbl);
 185         var += d2;
 186         skew += dbl * d2;
 187         kurtosis += d2 * d2;
 188     }
 189     var /= u.size();
 190     //double dev = std::sqrt(var);
 191     // In this case:
 192     //   skew     computes to 0./0. == nan
 193     //   kurtosis computes to 0./0. == nan
 194     //   x_skew     == inf
 195     //   x_kurtosis == inf
 196     //   These tests are commented out because UBSan warns about division by 0
 197 //    skew /= u.size() * dev * var;
 198 //    kurtosis /= u.size() * var * var;
 199 //    kurtosis -= 3;
 200     double x_mean = d.t() * d.p();
 201     double x_var = x_mean*(1-d.p());
 202 //    double x_skew = (1-2*d.p()) / std::sqrt(x_var);
 203 //    double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
 204     assert(mean == x_mean);
 205     assert(var == x_var);
 206 //    assert(skew == x_skew);
 207 //    assert(kurtosis == x_kurtosis);
 208 }
 209
 210 void
 211 test5()
 212 {
 213     typedef std::binomial_distribution<> D;
 214     typedef std::mt19937 G;
 215     G g;
 216     D d(40, 1);
 217     const int N = 100000;
 218     std::vector<D::result_type> u;
 219     for (int i = 0; i < N; ++i)
 220     {
 221         D::result_type v = d(g);
 222         assert(d.min() <= v && v <= d.max());
 223         u.push_back(v);
 224     }
 225     double mean = std::accumulate(u.begin(), u.end(),
 226                                           double(0)) / u.size();
 227     double var = 0;
 228     double skew = 0;
 229     double kurtosis = 0;
 230     for (unsigned i = 0; i < u.size(); ++i)
 231     {
 232         double dbl = (u[i] - mean);
 233         double d2 = sqr(dbl);
 234         var += d2;
 235         skew += dbl * d2;
 236         kurtosis += d2 * d2;
 237     }
 238     var /= u.size();
 239 //    double dev = std::sqrt(var);
 240     // In this case:
 241     //   skew     computes to 0./0. == nan
 242     //   kurtosis computes to 0./0. == nan
 243     //   x_skew     == -inf
 244     //   x_kurtosis == inf
 245     //   These tests are commented out because UBSan warns about division by 0
 246 //    skew /= u.size() * dev * var;
 247 //    kurtosis /= u.size() * var * var;
 248 //    kurtosis -= 3;
 249     double x_mean = d.t() * d.p();
 250     double x_var = x_mean*(1-d.p());
 251 //    double x_skew = (1-2*d.p()) / std::sqrt(x_var);
 252 //    double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
 253     assert(mean == x_mean);
 254     assert(var == x_var);
 255 //    assert(skew == x_skew);
 256 //    assert(kurtosis == x_kurtosis);
 257 }
 258
 259 void
 260 test6()
 261 {
 262     typedef std::binomial_distribution<> D;
 263     typedef std::mt19937 G;
 264     G g;
 265     D d(400, 0.5);
 266     const int N = 100000;
 267     std::vector<D::result_type> u;
 268     for (int i = 0; i < N; ++i)
 269     {
 270         D::result_type v = d(g);
 271         assert(d.min() <= v && v <= d.max());
 272         u.push_back(v);
 273     }
 274     double mean = std::accumulate(u.begin(), u.end(),
 275                                           double(0)) / u.size();
 276     double var = 0;
 277     double skew = 0;
 278     double kurtosis = 0;
 279     for (unsigned i = 0; i < u.size(); ++i)
 280     {
 281         double dbl = (u[i] - mean);
 282         double d2 = sqr(dbl);
 283         var += d2;
 284         skew += dbl * d2;
 285         kurtosis += d2 * d2;
 286     }
 287     var /= u.size();
 288     double dev = std::sqrt(var);
 289     skew /= u.size() * dev * var;
 290     kurtosis /= u.size() * var * var;
 291     kurtosis -= 3;
 292     double x_mean = d.t() * d.p();
 293     double x_var = x_mean*(1-d.p());
 294     double x_skew = (1-2*d.p()) / std::sqrt(x_var);
 295     double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
 296     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 297     assert(std::abs((var - x_var) / x_var) < 0.01);
 298     assert(std::abs(skew - x_skew) < 0.01);
 299     assert(std::abs(kurtosis - x_kurtosis) < 0.01);
 300 }
 301
 302 void
 303 test7()
 304 {
 305     typedef std::binomial_distribution<> D;
 306     typedef std::mt19937 G;
 307     G g;
 308     D d(1, 0.5);
 309     const int N = 100000;
 310     std::vector<D::result_type> u;
 311     for (int i = 0; i < N; ++i)
 312     {
 313         D::result_type v = d(g);
 314         assert(d.min() <= v && v <= d.max());
 315         u.push_back(v);
 316     }
 317     double mean = std::accumulate(u.begin(), u.end(),
 318                                           double(0)) / u.size();
 319     double var = 0;
 320     double skew = 0;
 321     double kurtosis = 0;
 322     for (unsigned i = 0; i < u.size(); ++i)
 323     {
 324         double dbl = (u[i] - mean);
 325         double d2 = sqr(dbl);
 326         var += d2;
 327         skew += dbl * d2;
 328         kurtosis += d2 * d2;
 329     }
 330     var /= u.size();
 331     double dev = std::sqrt(var);
 332     skew /= u.size() * dev * var;
 333     kurtosis /= u.size() * var * var;
 334     kurtosis -= 3;
 335     double x_mean = d.t() * d.p();
 336     double x_var = x_mean*(1-d.p());
 337     double x_skew = (1-2*d.p()) / std::sqrt(x_var);
 338     double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
 339     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 340     assert(std::abs((var - x_var) / x_var) < 0.01);
 341     assert(std::abs(skew - x_skew) < 0.01);
 342     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
 343 }
 344
 345 void
 346 test8()
 347 {
 348     const int N = 100000;
 349     std::mt19937 gen1;
 350     std::mt19937 gen2;
 351
 352     std::binomial_distribution<>         dist1(5, 0.1);
 353     std::binomial_distribution<unsigned> dist2(5, 0.1);
 354
 355     for(int i = 0; i < N; ++i) {
 356         int r1 = dist1(gen1);
 357         unsigned r2 = dist2(gen2);
 358         assert(r1 >= 0);
 359         assert(static_cast<unsigned>(r1) == r2);
 360     }
 361 }
 362
 363 void
 364 test9()
 365 {
 366     typedef std::binomial_distribution<> D;
 367     typedef std::mt19937 G;
 368     G g;
 369     D d(0, 0.005);
 370     const int N = 100000;
 371     std::vector<D::result_type> u;
 372     for (int i = 0; i < N; ++i)
 373     {
 374         D::result_type v = d(g);
 375         assert(d.min() <= v && v <= d.max());
 376         u.push_back(v);
 377     }
 378     double mean = std::accumulate(u.begin(), u.end(),
 379                                           double(0)) / u.size();
 380     double var = 0;
 381     double skew = 0;
 382     double kurtosis = 0;
 383     for (unsigned i = 0; i < u.size(); ++i)
 384     {
 385         double dbl = (u[i] - mean);
 386         double d2 = sqr(dbl);
 387         var += d2;
 388         skew += dbl * d2;
 389         kurtosis += d2 * d2;
 390     }
 391     var /= u.size();
 392 //    double dev = std::sqrt(var);
 393     // In this case:
 394     //   skew     computes to 0./0. == nan
 395     //   kurtosis computes to 0./0. == nan
 396     //   x_skew     == inf
 397     //   x_kurtosis == inf
 398     //   These tests are commented out because UBSan warns about division by 0
 399 //    skew /= u.size() * dev * var;
 400 //    kurtosis /= u.size() * var * var;
 401 //    kurtosis -= 3;
 402     double x_mean = d.t() * d.p();
 403     double x_var = x_mean*(1-d.p());
 404 //    double x_skew = (1-2*d.p()) / std::sqrt(x_var);
 405 //    double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
 406     assert(mean == x_mean);
 407     assert(var == x_var);
 408 //    assert(skew == x_skew);
 409 //    assert(kurtosis == x_kurtosis);
 410 }
 411
 412 void
 413 test10()
 414 {
 415     typedef std::binomial_distribution<> D;
 416     typedef std::mt19937 G;
 417     G g;
 418     D d(0, 0);
 419     const int N = 100000;
 420     std::vector<D::result_type> u;
 421     for (int i = 0; i < N; ++i)
 422     {
 423         D::result_type v = d(g);
 424         assert(d.min() <= v && v <= d.max());
 425         u.push_back(v);
 426     }
 427     double mean = std::accumulate(u.begin(), u.end(),
 428                                           double(0)) / u.size();
 429     double var = 0;
 430     double skew = 0;
 431     double kurtosis = 0;
 432     for (unsigned i = 0; i < u.size(); ++i)
 433     {
 434         double dbl = (u[i] - mean);
 435         double d2 = sqr(dbl);
 436         var += d2;
 437         skew += dbl * d2;
 438         kurtosis += d2 * d2;
 439     }
 440     var /= u.size();
 441 //    double dev = std::sqrt(var);
 442     // In this case:
 443     //   skew     computes to 0./0. == nan
 444     //   kurtosis computes to 0./0. == nan
 445     //   x_skew     == inf
 446     //   x_kurtosis == inf
 447     //   These tests are commented out because UBSan warns about division by 0
 448 //    skew /= u.size() * dev * var;
 449 //    kurtosis /= u.size() * var * var;
 450 //    kurtosis -= 3;
 451     double x_mean = d.t() * d.p();
 452     double x_var = x_mean*(1-d.p());
 453 //    double x_skew = (1-2*d.p()) / std::sqrt(x_var);
 454 //    double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
 455     assert(mean == x_mean);
 456     assert(var == x_var);
 457 //    assert(skew == x_skew);
 458 //    assert(kurtosis == x_kurtosis);
 459 }
 460
 461 void
 462 test11()
 463 {
 464     typedef std::binomial_distribution<> D;
 465     typedef std::mt19937 G;
 466     G g;
 467     D d(0, 1);
 468     const int N = 100000;
 469     std::vector<D::result_type> u;
 470     for (int i = 0; i < N; ++i)
 471     {
 472         D::result_type v = d(g);
 473         assert(d.min() <= v && v <= d.max());
 474         u.push_back(v);
 475     }
 476     double mean = std::accumulate(u.begin(), u.end(),
 477                                           double(0)) / u.size();
 478     double var = 0;
 479     double skew = 0;
 480     double kurtosis = 0;
 481     for (unsigned i = 0; i < u.size(); ++i)
 482     {
 483         double dbl = (u[i] - mean);
 484         double d2 = sqr(dbl);
 485         var += d2;
 486         skew += dbl * d2;
 487         kurtosis += d2 * d2;
 488     }
 489     var /= u.size();
 490 //    double dev = std::sqrt(var);
 491     // In this case:
 492     //   skew     computes to 0./0. == nan
 493     //   kurtosis computes to 0./0. == nan
 494     //   x_skew     == -inf
 495     //   x_kurtosis == inf
 496     //   These tests are commented out because UBSan warns about division by 0
 497 //    skew /= u.size() * dev * var;
 498 //    kurtosis /= u.size() * var * var;
 499 //    kurtosis -= 3;
 500     double x_mean = d.t() * d.p();
 501     double x_var = x_mean*(1-d.p());
 502 //    double x_skew = (1-2*d.p()) / std::sqrt(x_var);
 503 //    double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
 504     assert(mean == x_mean);
 505     assert(var == x_var);
 506 //    assert(skew == x_skew);
 507 //    assert(kurtosis == x_kurtosis);
 508 }
 509
 510 int main()
 511 {
 512     test1();
 513     test2();
 514     test3();
 515     test4();
 516     test5();
 517     test6();
 518     test7();
 519     test8();
 520     test9();
 521     test10();
 522     test11();
 523 }