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, const param_type& parm); |
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 D::param_type P; |
38 | typedef std::mt19937 G; |
39 | G g; |
40 | D d(-0.5, 1); |
41 | P p(0.5, 2); |
42 | const int N = 1000000; |
43 | std::vector<D::result_type> u; |
44 | for (int i = 0; i < N; ++i) |
45 | { |
46 | D::result_type v = d(g, p); |
47 | u.push_back(x: v); |
48 | } |
49 | double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size(); |
50 | double var = 0; |
51 | double skew = 0; |
52 | double kurtosis = 0; |
53 | for (unsigned i = 0; i < u.size(); ++i) |
54 | { |
55 | double dbl = (u[i] - mean); |
56 | double d2 = sqr(dbl); |
57 | var += d2; |
58 | skew += dbl * d2; |
59 | kurtosis += d2 * d2; |
60 | } |
61 | var /= u.size(); |
62 | double dev = std::sqrt(x: var); |
63 | skew /= u.size() * dev * var; |
64 | kurtosis /= u.size() * var * var; |
65 | kurtosis -= 3; |
66 | double x_mean = p.a() + p.b() * 0.577215665; |
67 | double x_var = sqr(p.b()) * 1.644934067; |
68 | double x_skew = 1.139547; |
69 | double x_kurtosis = 12./5; |
70 | assert(std::abs((mean - x_mean) / x_mean) < 0.01); |
71 | assert(std::abs((var - x_var) / x_var) < 0.01); |
72 | assert(std::abs((skew - x_skew) / x_skew) < 0.01); |
73 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); |
74 | } |
75 | |
76 | void |
77 | test2() |
78 | { |
79 | typedef std::extreme_value_distribution<> D; |
80 | typedef D::param_type P; |
81 | typedef std::mt19937 G; |
82 | G g; |
83 | D d(-0.5, 1); |
84 | P p(1, 2); |
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 | u.push_back(x: v); |
91 | } |
92 | double mean = std::accumulate(first: u.begin(), last: u.end(), init: 0.0) / u.size(); |
93 | double var = 0; |
94 | double skew = 0; |
95 | double kurtosis = 0; |
96 | for (unsigned i = 0; i < u.size(); ++i) |
97 | { |
98 | double dbl = (u[i] - mean); |
99 | double d2 = sqr(x: dbl); |
100 | var += d2; |
101 | skew += dbl * d2; |
102 | kurtosis += d2 * d2; |
103 | } |
104 | var /= u.size(); |
105 | double dev = std::sqrt(x: var); |
106 | skew /= u.size() * dev * var; |
107 | kurtosis /= u.size() * var * var; |
108 | kurtosis -= 3; |
109 | double x_mean = p.a() + p.b() * 0.577215665; |
110 | double x_var = sqr(x: p.b()) * 1.644934067; |
111 | double x_skew = 1.139547; |
112 | double x_kurtosis = 12./5; |
113 | assert(std::abs((mean - x_mean) / x_mean) < 0.01); |
114 | assert(std::abs((var - x_var) / x_var) < 0.01); |
115 | assert(std::abs((skew - x_skew) / x_skew) < 0.01); |
116 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); |
117 | } |
118 | |
119 | void |
120 | test3() |
121 | { |
122 | typedef std::extreme_value_distribution<> D; |
123 | typedef D::param_type P; |
124 | typedef std::mt19937 G; |
125 | G g; |
126 | D d(-0.5, 1); |
127 | P p(1.5, 3); |
128 | const int N = 1000000; |
129 | std::vector<D::result_type> u; |
130 | for (int i = 0; i < N; ++i) |
131 | { |
132 | D::result_type v = d(g, p); |
133 | u.push_back(x: v); |
134 | } |
135 | double mean = std::accumulate(first: u.begin(), last: u.end(), init: 0.0) / u.size(); |
136 | double var = 0; |
137 | double skew = 0; |
138 | double kurtosis = 0; |
139 | for (unsigned i = 0; i < u.size(); ++i) |
140 | { |
141 | double dbl = (u[i] - mean); |
142 | double d2 = sqr(x: dbl); |
143 | var += d2; |
144 | skew += dbl * d2; |
145 | kurtosis += d2 * d2; |
146 | } |
147 | var /= u.size(); |
148 | double dev = std::sqrt(x: var); |
149 | skew /= u.size() * dev * var; |
150 | kurtosis /= u.size() * var * var; |
151 | kurtosis -= 3; |
152 | double x_mean = p.a() + p.b() * 0.577215665; |
153 | double x_var = sqr(x: p.b()) * 1.644934067; |
154 | double x_skew = 1.139547; |
155 | double x_kurtosis = 12./5; |
156 | assert(std::abs((mean - x_mean) / x_mean) < 0.01); |
157 | assert(std::abs((var - x_var) / x_var) < 0.01); |
158 | assert(std::abs((skew - x_skew) / x_skew) < 0.01); |
159 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); |
160 | } |
161 | |
162 | void |
163 | test4() |
164 | { |
165 | typedef std::extreme_value_distribution<> D; |
166 | typedef D::param_type P; |
167 | typedef std::mt19937 G; |
168 | G g; |
169 | D d(-0.5, 1); |
170 | P p(3, 4); |
171 | const int N = 1000000; |
172 | std::vector<D::result_type> u; |
173 | for (int i = 0; i < N; ++i) |
174 | { |
175 | D::result_type v = d(g, p); |
176 | u.push_back(x: v); |
177 | } |
178 | double mean = std::accumulate(first: u.begin(), last: u.end(), init: 0.0) / u.size(); |
179 | double var = 0; |
180 | double skew = 0; |
181 | double kurtosis = 0; |
182 | for (unsigned i = 0; i < u.size(); ++i) |
183 | { |
184 | double dbl = (u[i] - mean); |
185 | double d2 = sqr(x: dbl); |
186 | var += d2; |
187 | skew += dbl * d2; |
188 | kurtosis += d2 * d2; |
189 | } |
190 | var /= u.size(); |
191 | double dev = std::sqrt(x: var); |
192 | skew /= u.size() * dev * var; |
193 | kurtosis /= u.size() * var * var; |
194 | kurtosis -= 3; |
195 | double x_mean = p.a() + p.b() * 0.577215665; |
196 | double x_var = sqr(x: p.b()) * 1.644934067; |
197 | double x_skew = 1.139547; |
198 | double x_kurtosis = 12./5; |
199 | assert(std::abs((mean - x_mean) / x_mean) < 0.01); |
200 | assert(std::abs((var - x_var) / x_var) < 0.01); |
201 | assert(std::abs((skew - x_skew) / x_skew) < 0.01); |
202 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); |
203 | } |
204 | |
205 | int main(int, char**) |
206 | { |
207 | test1(); |
208 | test2(); |
209 | test3(); |
210 | test4(); |
211 | |
212 | return 0; |
213 | } |
214 | |