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