xref: /llvm-project/libcxx/test/std/numerics/rand/rand.dist/rand.dist.bern/rand.dist.bern.bin/eval.pass.cpp (revision 09e3a360581dc36d0820d3fb6da9bd7cfed87b5d) 
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 <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 
28 template <class T>
29 T sqr(T x) {
30     return x * x;
31 }
32 
33 template <class T>
34 void 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 
75 template <class T>
76 void 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 
117 template <class T>
118 void 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 
159 template <class T>
160 void 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. == nan
190     //   kurtosis computes to 0./0. == nan
191     //   x_skew     == inf
192     //   x_kurtosis == inf
193     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 
208 template <class T>
209 void 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. == nan
239     //   kurtosis computes to 0./0. == nan
240     //   x_skew     == -inf
241     //   x_kurtosis == inf
242     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 
257 template <class T>
258 void 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 
299 template <class T>
300 void 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 
341 template <class T>
342 void 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 
359 template <class T>
360 void 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. == nan
390     //   kurtosis computes to 0./0. == nan
391     //   x_skew     == inf
392     //   x_kurtosis == inf
393     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 
408 template <class T>
409 void 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. == nan
439     //   kurtosis computes to 0./0. == nan
440     //   x_skew     == inf
441     //   x_kurtosis == inf
442     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 
457 template <class T>
458 void 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. == nan
488     //   kurtosis computes to 0./0. == nan
489     //   x_skew     == -inf
490     //   x_kurtosis == inf
491     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 
506 template <class T>
507 void 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 
521 int 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) // extension
533     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 #endif
539 #endif
540 
541     return 0;
542 }
543