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