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
24template <class T>
25inline
26T
27sqr(T x)
28{
29    return x * x;
30}
31
32void
33test1()
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
75void
76test2()
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
118void
119test3()
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
161void
162test4()
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
210void
211test5()
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
259void
260test6()
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
302void
303test7()
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
345void
346test8()
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
363void
364test9()
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
412void
413test10()
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
461void
462test11()
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
510int 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