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, const param_type& parm);
17
18#include <random>
19#include <numeric>
20#include <vector>
21#include <cassert>
22
23#include "test_macros.h"
24
25template <class T>
26inline
27T
28sqr(T x)
29{
30 return x * x;
31}
32
33int main(int, char**)
34{
35 {
36 typedef std::binomial_distribution<> D;
37 typedef D::param_type P;
38 typedef std::mt19937_64 G;
39 G g;
40 D d(16, .75);
41 P p(5, .75);
42 const int N = 1000000;
43 std::vector<D::result_type> u;
44 for (int i = 0; i < N; ++i)
45 {
46 D::result_type v = d(g, p);
47 assert(0 <= v && v <= p.t());
48 u.push_back(v);
49 }
50 double mean = std::accumulate(u.begin(), u.end(),
51 double(0)) / u.size();
52 double var = 0;
53 double skew = 0;
54 double kurtosis = 0;
55 for (unsigned i = 0; i < u.size(); ++i)
56 {
57 double dbl = (u[i] - mean);
58 double d2 = sqr(dbl);
59 var += d2;
60 skew += dbl * d2;
61 kurtosis += d2 * d2;
62 }
63 var /= u.size();
64 double dev = std::sqrt(x: var);
65 skew /= u.size() * dev * var;
66 kurtosis /= u.size() * var * var;
67 kurtosis -= 3;
68 double x_mean = p.t() * p.p();
69 double x_var = x_mean*(1-p.p());
70 double x_skew = (1-2*p.p()) / std::sqrt(x: x_var);
71 double x_kurtosis = (1-6*p.p()*(1-p.p())) / x_var;
72 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
73 assert(std::abs((var - x_var) / x_var) < 0.01);
74 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
75 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.04);
76 }
77 {
78 typedef std::binomial_distribution<> D;
79 typedef D::param_type P;
80 typedef std::mt19937 G;
81 G g;
82 D d(16, .75);
83 P p(30, .03125);
84 const int N = 100000;
85 std::vector<D::result_type> u;
86 for (int i = 0; i < N; ++i)
87 {
88 D::result_type v = d(g, p);
89 assert(0 <= v && v <= p.t());
90 u.push_back(v);
91 }
92 double mean = std::accumulate(u.begin(), u.end(),
93 double(0)) / u.size();
94 double var = 0;
95 double skew = 0;
96 double kurtosis = 0;
97 for (unsigned i = 0; i < u.size(); ++i)
98 {
99 double dbl = (u[i] - mean);
100 double d2 = sqr(dbl);
101 var += d2;
102 skew += dbl * d2;
103 kurtosis += d2 * d2;
104 }
105 var /= u.size();
106 double dev = std::sqrt(x: var);
107 skew /= u.size() * dev * var;
108 kurtosis /= u.size() * var * var;
109 kurtosis -= 3;
110 double x_mean = p.t() * p.p();
111 double x_var = x_mean*(1-p.p());
112 double x_skew = (1-2*p.p()) / std::sqrt(x: x_var);
113 double x_kurtosis = (1-6*p.p()*(1-p.p())) / x_var;
114 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
115 assert(std::abs((var - x_var) / x_var) < 0.01);
116 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
117 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
118 }
119 {
120 typedef std::binomial_distribution<> D;
121 typedef D::param_type P;
122 typedef std::mt19937 G;
123 G g;
124 D d(16, .75);
125 P p(40, .25);
126 const int N = 1000000;
127 std::vector<D::result_type> u;
128 for (int i = 0; i < N; ++i)
129 {
130 D::result_type v = d(g, p);
131 assert(0 <= v && v <= p.t());
132 u.push_back(v);
133 }
134 double mean = std::accumulate(u.begin(), u.end(),
135 double(0)) / u.size();
136 double var = 0;
137 double skew = 0;
138 double kurtosis = 0;
139 for (unsigned i = 0; i < u.size(); ++i)
140 {
141 double dbl = (u[i] - mean);
142 double d2 = sqr(dbl);
143 var += d2;
144 skew += dbl * d2;
145 kurtosis += d2 * d2;
146 }
147 var /= u.size();
148 double dev = std::sqrt(x: var);
149 skew /= u.size() * dev * var;
150 kurtosis /= u.size() * var * var;
151 kurtosis -= 3;
152 double x_mean = p.t() * p.p();
153 double x_var = x_mean*(1-p.p());
154 double x_skew = (1-2*p.p()) / std::sqrt(x: x_var);
155 double x_kurtosis = (1-6*p.p()*(1-p.p())) / x_var;
156 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
157 assert(std::abs((var - x_var) / x_var) < 0.01);
158 assert(std::abs((skew - x_skew) / x_skew) < 0.04);
159 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.3);
160 }
161
162 return 0;
163}
164

source code of libcxx/test/std/numerics/rand/rand.dist/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp