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, const param_type& parm);
17 
18 #include <random>
19 #include <cassert>
20 #include <cmath>
21 #include <numeric>
22 #include <vector>
23 
24 #include "test_macros.h"
25 
26 template <class T>
27 inline
28 T
29 sqr(T x)
30 {
31     return x * x;
32 }
33 
34 int main(int, char**)
35 {
36     {
37         typedef std::negative_binomial_distribution<> D;
38         typedef D::param_type P;
39         typedef std::minstd_rand G;
40         G g;
41         D d(16, .75);
42         P p(5, .75);
43         const int N = 1000000;
44         std::vector<D::result_type> u;
45         for (int i = 0; i < N; ++i)
46         {
47             D::result_type v = d(g, p);
48             assert(d.min() <= v && v <= d.max());
49             u.push_back(v);
50         }
51         double mean = std::accumulate(u.begin(), u.end(),
52                                               double(0)) / u.size();
53         double var = 0;
54         double skew = 0;
55         double kurtosis = 0;
56         for (unsigned i = 0; i < u.size(); ++i)
57         {
58             double dbl = (u[i] - mean);
59             double d2 = sqr(dbl);
60             var += d2;
61             skew += dbl * d2;
62             kurtosis += d2 * d2;
63         }
64         var /= u.size();
65         double dev = std::sqrt(var);
66         skew /= u.size() * dev * var;
67         kurtosis /= u.size() * var * var;
68         kurtosis -= 3;
69         double x_mean = p.k() * (1 - p.p()) / p.p();
70         double x_var = x_mean / p.p();
71         double x_skew = (2 - p.p()) / std::sqrt(p.k() * (1 - p.p()));
72         double x_kurtosis = 6. / p.k() + sqr(p.p()) / (p.k() * (1 - p.p()));
73         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
74         assert(std::abs((var - x_var) / x_var) < 0.01);
75         assert(std::abs((skew - x_skew) / x_skew) < 0.01);
76         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
77     }
78     {
79         typedef std::negative_binomial_distribution<> D;
80         typedef D::param_type P;
81         typedef std::mt19937 G;
82         G g;
83         D d(16, .75);
84         P p(30, .03125);
85         const int N = 1000000;
86         std::vector<D::result_type> u;
87         for (int i = 0; i < N; ++i)
88         {
89             D::result_type v = d(g, p);
90             assert(d.min() <= v && v <= d.max());
91             u.push_back(v);
92         }
93         double mean = std::accumulate(u.begin(), u.end(),
94                                               double(0)) / u.size();
95         double var = 0;
96         double skew = 0;
97         double kurtosis = 0;
98         for (unsigned i = 0; i < u.size(); ++i)
99         {
100             double dbl = (u[i] - mean);
101             double d2 = sqr(dbl);
102             var += d2;
103             skew += dbl * d2;
104             kurtosis += d2 * d2;
105         }
106         var /= u.size();
107         double dev = std::sqrt(var);
108         skew /= u.size() * dev * var;
109         kurtosis /= u.size() * var * var;
110         kurtosis -= 3;
111         double x_mean = p.k() * (1 - p.p()) / p.p();
112         double x_var = x_mean / p.p();
113         double x_skew = (2 - p.p()) / std::sqrt(p.k() * (1 - p.p()));
114         double x_kurtosis = 6. / p.k() + sqr(p.p()) / (p.k() * (1 - p.p()));
115         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
116         assert(std::abs((var - x_var) / x_var) < 0.01);
117         assert(std::abs((skew - x_skew) / x_skew) < 0.02);
118         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.1);
119     }
120     {
121         typedef std::negative_binomial_distribution<> D;
122         typedef D::param_type P;
123         typedef std::mt19937 G;
124         G g;
125         D d(16, .75);
126         P p(40, .25);
127         const int N = 1000000;
128         std::vector<D::result_type> u;
129         for (int i = 0; i < N; ++i)
130         {
131             D::result_type v = d(g, p);
132             assert(d.min() <= v && v <= d.max());
133             u.push_back(v);
134         }
135         double mean = std::accumulate(u.begin(), u.end(),
136                                               double(0)) / u.size();
137         double var = 0;
138         double skew = 0;
139         double kurtosis = 0;
140         for (unsigned i = 0; i < u.size(); ++i)
141         {
142             double dbl = (u[i] - mean);
143             double d2 = sqr(dbl);
144             var += d2;
145             skew += dbl * d2;
146             kurtosis += d2 * d2;
147         }
148         var /= u.size();
149         double dev = std::sqrt(var);
150         skew /= u.size() * dev * var;
151         kurtosis /= u.size() * var * var;
152         kurtosis -= 3;
153         double x_mean = p.k() * (1 - p.p()) / p.p();
154         double x_var = x_mean / p.p();
155         double x_skew = (2 - p.p()) / std::sqrt(p.k() * (1 - p.p()));
156         double x_kurtosis = 6. / p.k() + sqr(p.p()) / (p.k() * (1 - p.p()));
157         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
158         assert(std::abs((var - x_var) / x_var) < 0.01);
159         assert(std::abs((skew - x_skew) / x_skew) < 0.02);
160         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.08);
161     }
162 
163   return 0;
164 }
165