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