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
sqr(T x)27 sqr(T x)
28 {
29 return x * x;
30 }
31
32 void
test1()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
test2()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
test3()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
test4()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
test5()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
test6()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
test7()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
test8()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
test9()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
test10()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
test11()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
main()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