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
sqr(T x)28 sqr(T x)
29 {
30     return x * x;
31 }
32 
33 void
test1()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(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
test2()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(v);
95     }
96     double mean = std::accumulate(u.begin(), u.end(), 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(dbl);
104         var += d2;
105         skew += dbl * d2;
106         kurtosis += d2 * d2;
107     }
108     var /= u.size();
109     double dev = std::sqrt(var);
110     skew /= u.size() * dev * var;
111     kurtosis /= u.size() * var * var;
112     kurtosis -= 3;
113     double x_mean = std::exp(p.m() + sqr(p.s())/2);
114     double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
115     double x_skew = (std::exp(sqr(p.s())) + 2) *
116           std::sqrt((std::exp(sqr(p.s())) - 1));
117     double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
118                       3*std::exp(2*sqr(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
test3()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(v);
141     }
142     double mean = std::accumulate(u.begin(), u.end(), 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(dbl);
150         var += d2;
151         skew += dbl * d2;
152         kurtosis += d2 * d2;
153     }
154     var /= u.size();
155     double dev = std::sqrt(var);
156     skew /= u.size() * dev * var;
157     kurtosis /= u.size() * var * var;
158     kurtosis -= 3;
159     double x_mean = std::exp(p.m() + sqr(p.s())/2);
160     double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
161     double x_skew = (std::exp(sqr(p.s())) + 2) *
162           std::sqrt((std::exp(sqr(p.s())) - 1));
163     double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
164                       3*std::exp(2*sqr(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
test4()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(v);
187     }
188     double mean = std::accumulate(u.begin(), u.end(), 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(dbl);
196         var += d2;
197         skew += dbl * d2;
198         kurtosis += d2 * d2;
199     }
200     var /= u.size();
201     double dev = std::sqrt(var);
202     skew /= u.size() * dev * var;
203     kurtosis /= u.size() * var * var;
204     kurtosis -= 3;
205     double x_mean = std::exp(p.m() + sqr(p.s())/2);
206     double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
207     double x_skew = (std::exp(sqr(p.s())) + 2) *
208           std::sqrt((std::exp(sqr(p.s())) - 1));
209     double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
210                       3*std::exp(2*sqr(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
test5()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(v);
233     }
234     double mean = std::accumulate(u.begin(), u.end(), 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(dbl);
242         var += d2;
243         skew += dbl * d2;
244         kurtosis += d2 * d2;
245     }
246     var /= u.size();
247     double dev = std::sqrt(var);
248     skew /= u.size() * dev * var;
249     kurtosis /= u.size() * var * var;
250     kurtosis -= 3;
251     double x_mean = std::exp(p.m() + sqr(p.s())/2);
252     double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
253     double x_skew = (std::exp(sqr(p.s())) + 2) *
254           std::sqrt((std::exp(sqr(p.s())) - 1));
255     double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
256                       3*std::exp(2*sqr(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 
main(int,char **)263 int main(int, char**)
264 {
265     test1();
266     test2();
267     test3();
268     test4();
269     test5();
270 
271   return 0;
272 }
273