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