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 } 524
