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