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 negative_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 #include "test_macros.h"
24
25 template <class T>
26 inline
27 T
sqr(T x)28 sqr(T x)
29 {
30 return x * x;
31 }
32
33 void
test1()34 test1()
35 {
36 typedef std::negative_binomial_distribution<> D;
37 typedef std::minstd_rand G;
38 G g;
39 D d(5, .25);
40 const int N = 1000000;
41 std::vector<D::result_type> u;
42 for (int i = 0; i < N; ++i)
43 {
44 D::result_type v = d(g);
45 assert(d.min() <= v && v <= d.max());
46 u.push_back(v);
47 }
48 double mean = std::accumulate(u.begin(), u.end(),
49 double(0)) / u.size();
50 double var = 0;
51 double skew = 0;
52 double kurtosis = 0;
53 for (unsigned i = 0; i < u.size(); ++i)
54 {
55 double dbl = (u[i] - mean);
56 double d2 = sqr(dbl);
57 var += d2;
58 skew += dbl * d2;
59 kurtosis += d2 * d2;
60 }
61 var /= u.size();
62 double dev = std::sqrt(var);
63 skew /= u.size() * dev * var;
64 kurtosis /= u.size() * var * var;
65 kurtosis -= 3;
66 double x_mean = d.k() * (1 - d.p()) / d.p();
67 double x_var = x_mean / d.p();
68 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
69 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
70 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
71 assert(std::abs((var - x_var) / x_var) < 0.01);
72 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
73 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
74 }
75
76 void
test2()77 test2()
78 {
79 typedef std::negative_binomial_distribution<> D;
80 typedef std::mt19937 G;
81 G g;
82 D d(30, .03125);
83 const int N = 1000000;
84 std::vector<D::result_type> u;
85 for (int i = 0; i < N; ++i)
86 {
87 D::result_type v = d(g);
88 assert(d.min() <= v && v <= d.max());
89 u.push_back(v);
90 }
91 double mean = std::accumulate(u.begin(), u.end(),
92 double(0)) / u.size();
93 double var = 0;
94 double skew = 0;
95 double kurtosis = 0;
96 for (unsigned i = 0; i < u.size(); ++i)
97 {
98 double dbl = (u[i] - mean);
99 double d2 = sqr(dbl);
100 var += d2;
101 skew += dbl * d2;
102 kurtosis += d2 * d2;
103 }
104 var /= u.size();
105 double dev = std::sqrt(var);
106 skew /= u.size() * dev * var;
107 kurtosis /= u.size() * var * var;
108 kurtosis -= 3;
109 double x_mean = d.k() * (1 - d.p()) / d.p();
110 double x_var = x_mean / d.p();
111 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
112 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
113 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
114 assert(std::abs((var - x_var) / x_var) < 0.01);
115 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
116 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
117 }
118
119 void
test3()120 test3()
121 {
122 typedef std::negative_binomial_distribution<> D;
123 typedef std::mt19937 G;
124 G g;
125 D d(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);
131 assert(d.min() <= v && v <= d.max());
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(var);
149 skew /= u.size() * dev * var;
150 kurtosis /= u.size() * var * var;
151 kurtosis -= 3;
152 double x_mean = d.k() * (1 - d.p()) / d.p();
153 double x_var = x_mean / d.p();
154 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
155 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
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.01);
159 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
160 }
161
162 void
test4()163 test4()
164 {
165 typedef std::negative_binomial_distribution<> D;
166 typedef std::mt19937 G;
167 G g;
168 D d(40, 1);
169 const int N = 1000;
170 std::vector<D::result_type> u;
171 for (int i = 0; i < N; ++i)
172 {
173 D::result_type v = d(g);
174 assert(d.min() <= v && v <= d.max());
175 u.push_back(v);
176 }
177 double mean = std::accumulate(u.begin(), u.end(),
178 double(0)) / u.size();
179 double var = 0;
180 double skew = 0;
181 double kurtosis = 0;
182 for (unsigned i = 0; i < u.size(); ++i)
183 {
184 double dbl = (u[i] - mean);
185 double d2 = sqr(dbl);
186 var += d2;
187 skew += dbl * d2;
188 kurtosis += d2 * d2;
189 }
190 var /= u.size();
191 double dev = std::sqrt(var);
192 skew /= u.size() * dev * var;
193 kurtosis /= u.size() * var * var;
194 kurtosis -= 3;
195 double x_mean = d.k() * (1 - d.p()) / d.p();
196 double x_var = x_mean / d.p();
197 // double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
198 // double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
199 assert(mean == x_mean);
200 assert(var == x_var);
201 }
202
203 void
test5()204 test5()
205 {
206 typedef std::negative_binomial_distribution<> D;
207 typedef std::mt19937 G;
208 G g;
209 D d(400, 0.5);
210 const int N = 1000000;
211 std::vector<D::result_type> u;
212 for (int i = 0; i < N; ++i)
213 {
214 D::result_type v = d(g);
215 assert(d.min() <= v && v <= d.max());
216 u.push_back(v);
217 }
218 double mean = std::accumulate(u.begin(), u.end(),
219 double(0)) / u.size();
220 double var = 0;
221 double skew = 0;
222 double kurtosis = 0;
223 for (unsigned i = 0; i < u.size(); ++i)
224 {
225 double dbl = (u[i] - mean);
226 double d2 = sqr(dbl);
227 var += d2;
228 skew += dbl * d2;
229 kurtosis += d2 * d2;
230 }
231 var /= u.size();
232 double dev = std::sqrt(var);
233 skew /= u.size() * dev * var;
234 kurtosis /= u.size() * var * var;
235 kurtosis -= 3;
236 double x_mean = d.k() * (1 - d.p()) / d.p();
237 double x_var = x_mean / d.p();
238 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
239 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
240 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
241 assert(std::abs((var - x_var) / x_var) < 0.01);
242 assert(std::abs((skew - x_skew) / x_skew) < 0.04);
243 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.05);
244 }
245
246 void
test6()247 test6()
248 {
249 typedef std::negative_binomial_distribution<> D;
250 typedef std::mt19937 G;
251 G g;
252 D d(1, 0.05);
253 const int N = 1000000;
254 std::vector<D::result_type> u;
255 for (int i = 0; i < N; ++i)
256 {
257 D::result_type v = d(g);
258 assert(d.min() <= v && v <= d.max());
259 u.push_back(v);
260 }
261 double mean = std::accumulate(u.begin(), u.end(),
262 double(0)) / u.size();
263 double var = 0;
264 double skew = 0;
265 double kurtosis = 0;
266 for (unsigned i = 0; i < u.size(); ++i)
267 {
268 double dbl = (u[i] - mean);
269 double d2 = sqr(dbl);
270 var += d2;
271 skew += dbl * d2;
272 kurtosis += d2 * d2;
273 }
274 var /= u.size();
275 double dev = std::sqrt(var);
276 skew /= u.size() * dev * var;
277 kurtosis /= u.size() * var * var;
278 kurtosis -= 3;
279 double x_mean = d.k() * (1 - d.p()) / d.p();
280 double x_var = x_mean / d.p();
281 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
282 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
283 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
284 assert(std::abs((var - x_var) / x_var) < 0.01);
285 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
286 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
287 }
288
main(int,char **)289 int main(int, char**)
290 {
291 test1();
292 test2();
293 test3();
294 test4();
295 test5();
296 test6();
297
298 return 0;
299 }
300