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