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 lognormal_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::lognormal_distribution<> D; |
37 | typedef D::param_type P; |
38 | typedef std::mt19937 G; |
39 | G g; |
40 | D d; |
41 | P p(-1./8192, 0.015625); |
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 | assert(v > 0); |
48 | u.push_back(v); |
49 | } |
50 | double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size(); |
51 | double var = 0; |
52 | double skew = 0; |
53 | double kurtosis = 0; |
54 | for (unsigned i = 0; i < u.size(); ++i) |
55 | { |
56 | double dbl = (u[i] - mean); |
57 | double d2 = sqr(dbl); |
58 | var += d2; |
59 | skew += dbl * d2; |
60 | kurtosis += d2 * d2; |
61 | } |
62 | var /= u.size(); |
63 | double dev = std::sqrt(x: var); |
64 | skew /= u.size() * dev * var; |
65 | kurtosis /= u.size() * var * var; |
66 | kurtosis -= 3; |
67 | double x_mean = std::exp(p.m() + sqr(p.s())/2); |
68 | double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s())); |
69 | double x_skew = (std::exp(sqr(p.s())) + 2) * |
70 | std::sqrt((std::exp(sqr(p.s())) - 1)); |
71 | double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) + |
72 | 3*std::exp(2*sqr(p.s())) - 6; |
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.05); |
76 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.25); |
77 | } |
78 | |
79 | void |
80 | test2() |
81 | { |
82 | typedef std::lognormal_distribution<> D; |
83 | typedef D::param_type P; |
84 | typedef std::mt19937 G; |
85 | G g; |
86 | D d; |
87 | P p(-1./32, 0.25); |
88 | const int N = 1000000; |
89 | std::vector<D::result_type> u; |
90 | for (int i = 0; i < N; ++i) |
91 | { |
92 | D::result_type v = d(g, p); |
93 | assert(v > 0); |
94 | u.push_back(x: v); |
95 | } |
96 | double mean = std::accumulate(first: u.begin(), last: u.end(), init: 0.0) / u.size(); |
97 | double var = 0; |
98 | double skew = 0; |
99 | double kurtosis = 0; |
100 | for (unsigned i = 0; i < u.size(); ++i) |
101 | { |
102 | double dbl = (u[i] - mean); |
103 | double d2 = sqr(x: dbl); |
104 | var += d2; |
105 | skew += dbl * d2; |
106 | kurtosis += d2 * d2; |
107 | } |
108 | var /= u.size(); |
109 | double dev = std::sqrt(x: var); |
110 | skew /= u.size() * dev * var; |
111 | kurtosis /= u.size() * var * var; |
112 | kurtosis -= 3; |
113 | double x_mean = std::exp(x: p.m() + sqr(x: p.s())/2); |
114 | double x_var = (std::exp(x: sqr(x: p.s())) - 1) * std::exp(x: 2*p.m() + sqr(x: p.s())); |
115 | double x_skew = (std::exp(x: sqr(x: p.s())) + 2) * |
116 | std::sqrt(x: (std::exp(x: sqr(x: p.s())) - 1)); |
117 | double x_kurtosis = std::exp(x: 4*sqr(x: p.s())) + 2*std::exp(x: 3*sqr(x: p.s())) + |
118 | 3*std::exp(x: 2*sqr(x: p.s())) - 6; |
119 | assert(std::abs((mean - x_mean) / x_mean) < 0.01); |
120 | assert(std::abs((var - x_var) / x_var) < 0.01); |
121 | assert(std::abs((skew - x_skew) / x_skew) < 0.01); |
122 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03); |
123 | } |
124 | |
125 | void |
126 | test3() |
127 | { |
128 | typedef std::lognormal_distribution<> D; |
129 | typedef D::param_type P; |
130 | typedef std::mt19937 G; |
131 | G g; |
132 | D d; |
133 | P p(-1./8, 0.5); |
134 | const int N = 1000000; |
135 | std::vector<D::result_type> u; |
136 | for (int i = 0; i < N; ++i) |
137 | { |
138 | D::result_type v = d(g, p); |
139 | assert(v > 0); |
140 | u.push_back(x: v); |
141 | } |
142 | double mean = std::accumulate(first: u.begin(), last: u.end(), init: 0.0) / u.size(); |
143 | double var = 0; |
144 | double skew = 0; |
145 | double kurtosis = 0; |
146 | for (unsigned i = 0; i < u.size(); ++i) |
147 | { |
148 | double dbl = (u[i] - mean); |
149 | double d2 = sqr(x: dbl); |
150 | var += d2; |
151 | skew += dbl * d2; |
152 | kurtosis += d2 * d2; |
153 | } |
154 | var /= u.size(); |
155 | double dev = std::sqrt(x: var); |
156 | skew /= u.size() * dev * var; |
157 | kurtosis /= u.size() * var * var; |
158 | kurtosis -= 3; |
159 | double x_mean = std::exp(x: p.m() + sqr(x: p.s())/2); |
160 | double x_var = (std::exp(x: sqr(x: p.s())) - 1) * std::exp(x: 2*p.m() + sqr(x: p.s())); |
161 | double x_skew = (std::exp(x: sqr(x: p.s())) + 2) * |
162 | std::sqrt(x: (std::exp(x: sqr(x: p.s())) - 1)); |
163 | double x_kurtosis = std::exp(x: 4*sqr(x: p.s())) + 2*std::exp(x: 3*sqr(x: p.s())) + |
164 | 3*std::exp(x: 2*sqr(x: p.s())) - 6; |
165 | assert(std::abs((mean - x_mean) / x_mean) < 0.01); |
166 | assert(std::abs((var - x_var) / x_var) < 0.01); |
167 | assert(std::abs((skew - x_skew) / x_skew) < 0.02); |
168 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.05); |
169 | } |
170 | |
171 | void |
172 | test4() |
173 | { |
174 | typedef std::lognormal_distribution<> D; |
175 | typedef D::param_type P; |
176 | typedef std::mt19937 G; |
177 | G g; |
178 | D d(3, 4); |
179 | P p; |
180 | const int N = 1000000; |
181 | std::vector<D::result_type> u; |
182 | for (int i = 0; i < N; ++i) |
183 | { |
184 | D::result_type v = d(g, p); |
185 | assert(v > 0); |
186 | u.push_back(x: v); |
187 | } |
188 | double mean = std::accumulate(first: u.begin(), last: u.end(), init: 0.0) / u.size(); |
189 | double var = 0; |
190 | double skew = 0; |
191 | double kurtosis = 0; |
192 | for (unsigned i = 0; i < u.size(); ++i) |
193 | { |
194 | double dbl = (u[i] - mean); |
195 | double d2 = sqr(x: dbl); |
196 | var += d2; |
197 | skew += dbl * d2; |
198 | kurtosis += d2 * d2; |
199 | } |
200 | var /= u.size(); |
201 | double dev = std::sqrt(x: var); |
202 | skew /= u.size() * dev * var; |
203 | kurtosis /= u.size() * var * var; |
204 | kurtosis -= 3; |
205 | double x_mean = std::exp(x: p.m() + sqr(x: p.s())/2); |
206 | double x_var = (std::exp(x: sqr(x: p.s())) - 1) * std::exp(x: 2*p.m() + sqr(x: p.s())); |
207 | double x_skew = (std::exp(x: sqr(x: p.s())) + 2) * |
208 | std::sqrt(x: (std::exp(x: sqr(x: p.s())) - 1)); |
209 | double x_kurtosis = std::exp(x: 4*sqr(x: p.s())) + 2*std::exp(x: 3*sqr(x: p.s())) + |
210 | 3*std::exp(x: 2*sqr(x: p.s())) - 6; |
211 | assert(std::abs((mean - x_mean) / x_mean) < 0.01); |
212 | assert(std::abs((var - x_var) / x_var) < 0.02); |
213 | assert(std::abs((skew - x_skew) / x_skew) < 0.08); |
214 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.4); |
215 | } |
216 | |
217 | void |
218 | test5() |
219 | { |
220 | typedef std::lognormal_distribution<> D; |
221 | typedef D::param_type P; |
222 | typedef std::mt19937 G; |
223 | G g; |
224 | D d; |
225 | P p(-0.78125, 1.25); |
226 | const int N = 1000000; |
227 | std::vector<D::result_type> u; |
228 | for (int i = 0; i < N; ++i) |
229 | { |
230 | D::result_type v = d(g, p); |
231 | assert(v > 0); |
232 | u.push_back(x: v); |
233 | } |
234 | double mean = std::accumulate(first: u.begin(), last: u.end(), init: 0.0) / u.size(); |
235 | double var = 0; |
236 | double skew = 0; |
237 | double kurtosis = 0; |
238 | for (unsigned i = 0; i < u.size(); ++i) |
239 | { |
240 | double dbl = (u[i] - mean); |
241 | double d2 = sqr(x: dbl); |
242 | var += d2; |
243 | skew += dbl * d2; |
244 | kurtosis += d2 * d2; |
245 | } |
246 | var /= u.size(); |
247 | double dev = std::sqrt(x: var); |
248 | skew /= u.size() * dev * var; |
249 | kurtosis /= u.size() * var * var; |
250 | kurtosis -= 3; |
251 | double x_mean = std::exp(x: p.m() + sqr(x: p.s())/2); |
252 | double x_var = (std::exp(x: sqr(x: p.s())) - 1) * std::exp(x: 2*p.m() + sqr(x: p.s())); |
253 | double x_skew = (std::exp(x: sqr(x: p.s())) + 2) * |
254 | std::sqrt(x: (std::exp(x: sqr(x: p.s())) - 1)); |
255 | double x_kurtosis = std::exp(x: 4*sqr(x: p.s())) + 2*std::exp(x: 3*sqr(x: p.s())) + |
256 | 3*std::exp(x: 2*sqr(x: p.s())) - 6; |
257 | assert(std::abs((mean - x_mean) / x_mean) < 0.01); |
258 | assert(std::abs((var - x_var) / x_var) < 0.04); |
259 | assert(std::abs((skew - x_skew) / x_skew) < 0.2); |
260 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.7); |
261 | } |
262 | |
263 | int main(int, char**) |
264 | { |
265 | test1(); |
266 | test2(); |
267 | test3(); |
268 | test4(); |
269 | test5(); |
270 | |
271 | return 0; |
272 | } |
273 | |