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 negative_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
24 template <class T>
25 inline
26 T
sqr(T x)27 sqr(T x)
28 {
29 return x * x;
30 }
31
32 void
test1()33 test1()
34 {
35 typedef std::negative_binomial_distribution<> D;
36 typedef std::minstd_rand G;
37 G g;
38 D d(5, .25);
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.k() * (1 - d.p()) / d.p();
66 double x_var = x_mean / d.p();
67 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
68 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
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.02);
73 }
74
75 void
test2()76 test2()
77 {
78 typedef std::negative_binomial_distribution<> D;
79 typedef std::mt19937 G;
80 G g;
81 D d(30, .03125);
82 const int N = 1000000;
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.k() * (1 - d.p()) / d.p();
109 double x_var = x_mean / d.p();
110 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
111 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
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
118 void
test3()119 test3()
120 {
121 typedef std::negative_binomial_distribution<> D;
122 typedef std::mt19937 G;
123 G g;
124 D d(40, .25);
125 const int N = 1000000;
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.k() * (1 - d.p()) / d.p();
152 double x_var = x_mean / d.p();
153 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
154 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
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.01);
158 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
159 }
160
161 void
test4()162 test4()
163 {
164 typedef std::negative_binomial_distribution<> D;
165 typedef std::mt19937 G;
166 G g;
167 D d(40, 1);
168 const int N = 1000;
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 skew /= u.size() * dev * var;
192 kurtosis /= u.size() * var * var;
193 kurtosis -= 3;
194 double x_mean = d.k() * (1 - d.p()) / d.p();
195 double x_var = x_mean / d.p();
196 // double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
197 // double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
198 assert(mean == x_mean);
199 assert(var == x_var);
200 }
201
202 void
test5()203 test5()
204 {
205 typedef std::negative_binomial_distribution<> D;
206 typedef std::mt19937 G;
207 G g;
208 D d(400, 0.5);
209 const int N = 1000000;
210 std::vector<D::result_type> u;
211 for (int i = 0; i < N; ++i)
212 {
213 D::result_type v = d(g);
214 assert(d.min() <= v && v <= d.max());
215 u.push_back(v);
216 }
217 double mean = std::accumulate(u.begin(), u.end(),
218 double(0)) / u.size();
219 double var = 0;
220 double skew = 0;
221 double kurtosis = 0;
222 for (unsigned i = 0; i < u.size(); ++i)
223 {
224 double dbl = (u[i] - mean);
225 double d2 = sqr(dbl);
226 var += d2;
227 skew += dbl * d2;
228 kurtosis += d2 * d2;
229 }
230 var /= u.size();
231 double dev = std::sqrt(var);
232 skew /= u.size() * dev * var;
233 kurtosis /= u.size() * var * var;
234 kurtosis -= 3;
235 double x_mean = d.k() * (1 - d.p()) / d.p();
236 double x_var = x_mean / d.p();
237 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
238 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
239 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
240 assert(std::abs((var - x_var) / x_var) < 0.01);
241 assert(std::abs((skew - x_skew) / x_skew) < 0.04);
242 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.05);
243 }
244
245 void
test6()246 test6()
247 {
248 typedef std::negative_binomial_distribution<> D;
249 typedef std::mt19937 G;
250 G g;
251 D d(1, 0.05);
252 const int N = 1000000;
253 std::vector<D::result_type> u;
254 for (int i = 0; i < N; ++i)
255 {
256 D::result_type v = d(g);
257 assert(d.min() <= v && v <= d.max());
258 u.push_back(v);
259 }
260 double mean = std::accumulate(u.begin(), u.end(),
261 double(0)) / u.size();
262 double var = 0;
263 double skew = 0;
264 double kurtosis = 0;
265 for (unsigned i = 0; i < u.size(); ++i)
266 {
267 double dbl = (u[i] - mean);
268 double d2 = sqr(dbl);
269 var += d2;
270 skew += dbl * d2;
271 kurtosis += d2 * d2;
272 }
273 var /= u.size();
274 double dev = std::sqrt(var);
275 skew /= u.size() * dev * var;
276 kurtosis /= u.size() * var * var;
277 kurtosis -= 3;
278 double x_mean = d.k() * (1 - d.p()) / d.p();
279 double x_var = x_mean / d.p();
280 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
281 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
282 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
283 assert(std::abs((var - x_var) / x_var) < 0.01);
284 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
285 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
286 }
287
main()288 int main()
289 {
290 test1();
291 test2();
292 test3();
293 test4();
294 test5();
295 test6();
296 }
297