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
24template <class T>
25inline
26T
27sqr(T x)
28{
29    return x * x;
30}
31
32void
33test1()
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
75void
76test2()
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
118void
119test3()
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
161void
162test4()
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
202void
203test5()
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
245void
246test6()
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
288int main()
289{
290    test1();
291    test2();
292    test3();
293    test4();
294    test5();
295    test6();
296}
297