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 RealType = double> |
14 | // class extreme_value_distribution |
15 | |
16 | // template<class _URNG> result_type operator()(_URNG& g); |
17 | |
18 | #include <random> |
19 | #include <cassert> |
20 | #include <vector> |
21 | #include <numeric> |
22 | |
23 | #include "test_macros.h" |
24 | |
25 | template <class T> |
26 | inline |
27 | T |
28 | sqr(T x) |
29 | { |
30 | return x * x; |
31 | } |
32 | |
33 | void |
34 | test1() |
35 | { |
36 | typedef std::extreme_value_distribution<> D; |
37 | typedef std::mt19937 G; |
38 | G g; |
39 | D d(0.5, 2); |
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 | u.push_back(x: v); |
46 | } |
47 | double mean = std::accumulate(u.begin(), u.end(), 0.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(x: var); |
61 | skew /= u.size() * dev * var; |
62 | kurtosis /= u.size() * var * var; |
63 | kurtosis -= 3; |
64 | double x_mean = d.a() + d.b() * 0.577215665; |
65 | double x_var = sqr(d.b()) * 1.644934067; |
66 | double x_skew = 1.139547; |
67 | double x_kurtosis = 12./5; |
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.01); |
72 | } |
73 | |
74 | void |
75 | test2() |
76 | { |
77 | typedef std::extreme_value_distribution<> D; |
78 | typedef std::mt19937 G; |
79 | G g; |
80 | D d(1, 2); |
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 | u.push_back(x: v); |
87 | } |
88 | double mean = std::accumulate(first: u.begin(), last: u.end(), init: 0.0) / u.size(); |
89 | double var = 0; |
90 | double skew = 0; |
91 | double kurtosis = 0; |
92 | for (unsigned i = 0; i < u.size(); ++i) |
93 | { |
94 | double dbl = (u[i] - mean); |
95 | double d2 = sqr(x: dbl); |
96 | var += d2; |
97 | skew += dbl * d2; |
98 | kurtosis += d2 * d2; |
99 | } |
100 | var /= u.size(); |
101 | double dev = std::sqrt(x: var); |
102 | skew /= u.size() * dev * var; |
103 | kurtosis /= u.size() * var * var; |
104 | kurtosis -= 3; |
105 | double x_mean = d.a() + d.b() * 0.577215665; |
106 | double x_var = sqr(x: d.b()) * 1.644934067; |
107 | double x_skew = 1.139547; |
108 | double x_kurtosis = 12./5; |
109 | assert(std::abs((mean - x_mean) / x_mean) < 0.01); |
110 | assert(std::abs((var - x_var) / x_var) < 0.01); |
111 | assert(std::abs((skew - x_skew) / x_skew) < 0.01); |
112 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); |
113 | } |
114 | |
115 | void |
116 | test3() |
117 | { |
118 | typedef std::extreme_value_distribution<> D; |
119 | typedef std::mt19937 G; |
120 | G g; |
121 | D d(1.5, 3); |
122 | const int N = 1000000; |
123 | std::vector<D::result_type> u; |
124 | for (int i = 0; i < N; ++i) |
125 | { |
126 | D::result_type v = d(g); |
127 | u.push_back(x: v); |
128 | } |
129 | double mean = std::accumulate(first: u.begin(), last: u.end(), init: 0.0) / u.size(); |
130 | double var = 0; |
131 | double skew = 0; |
132 | double kurtosis = 0; |
133 | for (unsigned i = 0; i < u.size(); ++i) |
134 | { |
135 | double dbl = (u[i] - mean); |
136 | double d2 = sqr(x: dbl); |
137 | var += d2; |
138 | skew += dbl * d2; |
139 | kurtosis += d2 * d2; |
140 | } |
141 | var /= u.size(); |
142 | double dev = std::sqrt(x: var); |
143 | skew /= u.size() * dev * var; |
144 | kurtosis /= u.size() * var * var; |
145 | kurtosis -= 3; |
146 | double x_mean = d.a() + d.b() * 0.577215665; |
147 | double x_var = sqr(x: d.b()) * 1.644934067; |
148 | double x_skew = 1.139547; |
149 | double x_kurtosis = 12./5; |
150 | assert(std::abs((mean - x_mean) / x_mean) < 0.01); |
151 | assert(std::abs((var - x_var) / x_var) < 0.01); |
152 | assert(std::abs((skew - x_skew) / x_skew) < 0.01); |
153 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); |
154 | } |
155 | |
156 | void |
157 | test4() |
158 | { |
159 | typedef std::extreme_value_distribution<> D; |
160 | typedef std::mt19937 G; |
161 | G g; |
162 | D d(3, 4); |
163 | const int N = 1000000; |
164 | std::vector<D::result_type> u; |
165 | for (int i = 0; i < N; ++i) |
166 | { |
167 | D::result_type v = d(g); |
168 | u.push_back(x: v); |
169 | } |
170 | double mean = std::accumulate(first: u.begin(), last: u.end(), init: 0.0) / u.size(); |
171 | double var = 0; |
172 | double skew = 0; |
173 | double kurtosis = 0; |
174 | for (unsigned i = 0; i < u.size(); ++i) |
175 | { |
176 | double dbl = (u[i] - mean); |
177 | double d2 = sqr(x: dbl); |
178 | var += d2; |
179 | skew += dbl * d2; |
180 | kurtosis += d2 * d2; |
181 | } |
182 | var /= u.size(); |
183 | double dev = std::sqrt(x: var); |
184 | skew /= u.size() * dev * var; |
185 | kurtosis /= u.size() * var * var; |
186 | kurtosis -= 3; |
187 | double x_mean = d.a() + d.b() * 0.577215665; |
188 | double x_var = sqr(x: d.b()) * 1.644934067; |
189 | double x_skew = 1.139547; |
190 | double x_kurtosis = 12./5; |
191 | assert(std::abs((mean - x_mean) / x_mean) < 0.01); |
192 | assert(std::abs((var - x_var) / x_var) < 0.01); |
193 | assert(std::abs((skew - x_skew) / x_skew) < 0.01); |
194 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); |
195 | } |
196 | |
197 | int main(int, char**) |
198 | { |
199 | test1(); |
200 | test2(); |
201 | test3(); |
202 | test4(); |
203 | |
204 | return 0; |
205 | } |
206 | |