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// This test is super slow, in particular with msan or tsan. In order to avoid timeouts and to
12// spend less time waiting for this particular test to complete we compile with optimizations.
13// ADDITIONAL_COMPILE_FLAGS(msan): -O1
14// ADDITIONAL_COMPILE_FLAGS(tsan): -O1
15
16// FIXME: This and other tests fail under GCC with optimizations enabled.
17// More investigation is needed, but it appears that GCC is performing more constant folding.
18
19// <random>
20
21// template<class IntType = int>
22// class negative_binomial_distribution
23
24// template<class _URNG> result_type operator()(_URNG& g);
25
26#include <random>
27#include <cassert>
28#include <cmath>
29#include <numeric>
30#include <vector>
31
32#include "test_macros.h"
33
34template <class T>
35T sqr(T x) {
36 return x * x;
37}
38
39template <class T>
40void test1() {
41 typedef std::negative_binomial_distribution<T> D;
42 typedef std::minstd_rand G;
43 G g;
44 D d(5, .25);
45 const int N = 1000000;
46 std::vector<typename D::result_type> u;
47 for (int i = 0; i < N; ++i)
48 {
49 typename D::result_type v = d(g);
50 assert(d.min() <= v && v <= d.max());
51 u.push_back(v);
52 }
53 double mean = std::accumulate(u.begin(), u.end(),
54 double(0)) / u.size();
55 double var = 0;
56 double skew = 0;
57 double kurtosis = 0;
58 for (unsigned i = 0; i < u.size(); ++i)
59 {
60 double dbl = (u[i] - mean);
61 double d2 = sqr(x: dbl);
62 var += d2;
63 skew += dbl * d2;
64 kurtosis += d2 * d2;
65 }
66 var /= u.size();
67 double dev = std::sqrt(x: var);
68 skew /= u.size() * dev * var;
69 kurtosis /= u.size() * var * var;
70 kurtosis -= 3;
71 double x_mean = d.k() * (1 - d.p()) / d.p();
72 double x_var = x_mean / d.p();
73 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
74 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
75 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
76 assert(std::abs((var - x_var) / x_var) < 0.01);
77 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
78 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
79}
80
81template <class T>
82void test2() {
83 typedef std::negative_binomial_distribution<T> D;
84 typedef std::mt19937 G;
85 G g;
86 D d(30, .03125);
87 const int N = 1000000;
88 std::vector<typename D::result_type> u;
89 for (int i = 0; i < N; ++i)
90 {
91 typename D::result_type v = d(g);
92 assert(d.min() <= v && v <= d.max());
93 u.push_back(v);
94 }
95 double mean = std::accumulate(u.begin(), u.end(),
96 double(0)) / u.size();
97 double var = 0;
98 double skew = 0;
99 double kurtosis = 0;
100 for (unsigned i = 0; i < u.size(); ++i)
101 {
102 double dbl = (u[i] - mean);
103 double d2 = sqr(x: dbl);
104 var += d2;
105 skew += dbl * d2;
106 kurtosis += d2 * d2;
107 }
108 var /= u.size();
109 double dev = std::sqrt(x: var);
110 skew /= u.size() * dev * var;
111 kurtosis /= u.size() * var * var;
112 kurtosis -= 3;
113 double x_mean = d.k() * (1 - d.p()) / d.p();
114 double x_var = x_mean / d.p();
115 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
116 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
117 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
118 assert(std::abs((var - x_var) / x_var) < 0.01);
119 assert(std::abs((skew - x_skew) / x_skew) < 0.02);
120 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.1);
121}
122
123template <class T>
124void test3() {
125 typedef std::negative_binomial_distribution<T> D;
126 typedef std::mt19937 G;
127 G g;
128 D d(40, .25);
129 const int N = 1000000;
130 std::vector<typename D::result_type> u;
131 for (int i = 0; i < N; ++i)
132 {
133 typename D::result_type v = d(g);
134 assert(d.min() <= v && v <= d.max());
135 u.push_back(v);
136 }
137 double mean = std::accumulate(u.begin(), u.end(),
138 double(0)) / u.size();
139 double var = 0;
140 double skew = 0;
141 double kurtosis = 0;
142 for (unsigned i = 0; i < u.size(); ++i)
143 {
144 double dbl = (u[i] - mean);
145 double d2 = sqr(x: dbl);
146 var += d2;
147 skew += dbl * d2;
148 kurtosis += d2 * d2;
149 }
150 var /= u.size();
151 double dev = std::sqrt(x: var);
152 skew /= u.size() * dev * var;
153 kurtosis /= u.size() * var * var;
154 kurtosis -= 3;
155 double x_mean = d.k() * (1 - d.p()) / d.p();
156 double x_var = x_mean / d.p();
157 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
158 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
159 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
160 assert(std::abs((var - x_var) / x_var) < 0.01);
161 assert(std::abs((skew - x_skew) / x_skew) < 0.02);
162 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.08);
163}
164
165template <class T>
166void test4() {
167 typedef std::negative_binomial_distribution<T> D;
168 typedef std::mt19937 G;
169 G g;
170 D d(40, 1);
171 const int N = 1000;
172 std::vector<typename D::result_type> u;
173 for (int i = 0; i < N; ++i)
174 {
175 typename D::result_type v = d(g);
176 assert(d.min() <= v && v <= d.max());
177 u.push_back(v);
178 }
179 double mean = std::accumulate(u.begin(), u.end(),
180 double(0)) / u.size();
181 double var = 0;
182 double skew = 0;
183 double kurtosis = 0;
184 for (unsigned i = 0; i < u.size(); ++i)
185 {
186 double dbl = (u[i] - mean);
187 double d2 = sqr(x: dbl);
188 var += d2;
189 skew += dbl * d2;
190 kurtosis += d2 * d2;
191 }
192 var /= u.size();
193 double dev = std::sqrt(x: var);
194 skew /= u.size() * dev * var;
195 kurtosis /= u.size() * var * var;
196 kurtosis -= 3;
197 double x_mean = d.k() * (1 - d.p()) / d.p();
198 double x_var = x_mean / d.p();
199 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
200 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
201 assert(mean == x_mean);
202 assert(var == x_var);
203 // assert(skew == x_skew);
204 (void)skew; (void)x_skew;
205 // assert(kurtosis == x_kurtosis);
206 (void)kurtosis; (void)x_kurtosis;
207}
208
209template <class T>
210void test5() {
211 typedef std::negative_binomial_distribution<T> D;
212 typedef std::mt19937 G;
213 G g;
214 D d(127, 0.5);
215 const int N = 1000000;
216 std::vector<typename D::result_type> u;
217 for (int i = 0; i < N; ++i)
218 {
219 typename D::result_type v = d(g);
220 assert(d.min() <= v && v <= d.max());
221 u.push_back(v);
222 }
223 double mean = std::accumulate(u.begin(), u.end(),
224 double(0)) / u.size();
225 double var = 0;
226 double skew = 0;
227 double kurtosis = 0;
228 for (unsigned i = 0; i < u.size(); ++i)
229 {
230 double dbl = (u[i] - mean);
231 double d2 = sqr(x: dbl);
232 var += d2;
233 skew += dbl * d2;
234 kurtosis += d2 * d2;
235 }
236 var /= u.size();
237 double dev = std::sqrt(x: var);
238 skew /= u.size() * dev * var;
239 kurtosis /= u.size() * var * var;
240 kurtosis -= 3;
241 double x_mean = d.k() * (1 - d.p()) / d.p();
242 double x_var = x_mean / d.p();
243 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
244 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
245 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
246 assert(std::abs((var - x_var) / x_var) < 0.01);
247 assert(std::abs((skew - x_skew) / x_skew) < 0.02);
248 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.3);
249}
250
251template <class T>
252void test6() {
253 typedef std::negative_binomial_distribution<T> D;
254 typedef std::mt19937 G;
255 G g;
256 D d(1, 0.05);
257 const int N = 1000000;
258 std::vector<typename D::result_type> u;
259 for (int i = 0; i < N; ++i)
260 {
261 typename D::result_type v = d(g);
262 assert(d.min() <= v && v <= d.max());
263 u.push_back(v);
264 }
265 double mean = std::accumulate(u.begin(), u.end(),
266 double(0)) / u.size();
267 double var = 0;
268 double skew = 0;
269 double kurtosis = 0;
270 for (unsigned i = 0; i < u.size(); ++i)
271 {
272 double dbl = (u[i] - mean);
273 double d2 = sqr(x: dbl);
274 var += d2;
275 skew += dbl * d2;
276 kurtosis += d2 * d2;
277 }
278 var /= u.size();
279 double dev = std::sqrt(x: var);
280 skew /= u.size() * dev * var;
281 kurtosis /= u.size() * var * var;
282 kurtosis -= 3;
283 double x_mean = d.k() * (1 - d.p()) / d.p();
284 double x_var = x_mean / d.p();
285 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
286 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
287 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
288 assert(std::abs((var - x_var) / x_var) < 0.01);
289 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
290 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
291}
292
293template <class T>
294void tests() {
295 test1<T>();
296 test2<T>();
297 test3<T>();
298 test4<T>();
299 test5<T>();
300 test6<T>();
301}
302
303int main(int, char**) {
304 tests<short>();
305 tests<int>();
306 tests<long>();
307 tests<long long>();
308
309 tests<unsigned short>();
310 tests<unsigned int>();
311 tests<unsigned long>();
312 tests<unsigned long long>();
313
314#if defined(_LIBCPP_VERSION) // extension
315 // TODO: std::negative_binomial_distribution currently doesn't work reliably with small types.
316 // tests<int8_t>();
317 // tests<uint8_t>();
318#if !defined(TEST_HAS_NO_INT128)
319 tests<__int128_t>();
320 tests<__uint128_t>();
321#endif
322#endif
323
324 return 0;
325}
326

source code of libcxx/test/std/numerics/rand/rand.dist/rand.dist.bern/rand.dist.bern.negbin/eval.pass.cpp