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
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 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(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
test2()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(v);
91     }
92     double mean = std::accumulate(u.begin(), u.end(), 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(dbl);
100         var += d2;
101         skew += dbl * d2;
102         kurtosis += d2 * d2;
103     }
104     var /= u.size();
105     double dev = std::sqrt(var);
106     skew /= u.size() * dev * var;
107     kurtosis /= u.size() * var * var;
108     kurtosis -= 3;
109     double x_mean = std::exp(d.m() + sqr(d.s())/2);
110     double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
111     double x_skew = (std::exp(sqr(d.s())) + 2) *
112           std::sqrt((std::exp(sqr(d.s())) - 1));
113     double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
114                       3*std::exp(2*sqr(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
test3()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(v);
135     }
136     double mean = std::accumulate(u.begin(), u.end(), 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(dbl);
144         var += d2;
145         skew += dbl * d2;
146         kurtosis += d2 * d2;
147     }
148     var /= u.size();
149     double dev = std::sqrt(var);
150     skew /= u.size() * dev * var;
151     kurtosis /= u.size() * var * var;
152     kurtosis -= 3;
153     double x_mean = std::exp(d.m() + sqr(d.s())/2);
154     double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
155     double x_skew = (std::exp(sqr(d.s())) + 2) *
156           std::sqrt((std::exp(sqr(d.s())) - 1));
157     double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
158                       3*std::exp(2*sqr(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
test4()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(v);
179     }
180     double mean = std::accumulate(u.begin(), u.end(), 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(dbl);
188         var += d2;
189         skew += dbl * d2;
190         kurtosis += d2 * d2;
191     }
192     var /= u.size();
193     double dev = std::sqrt(var);
194     skew /= u.size() * dev * var;
195     kurtosis /= u.size() * var * var;
196     kurtosis -= 3;
197     double x_mean = std::exp(d.m() + sqr(d.s())/2);
198     double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
199     double x_skew = (std::exp(sqr(d.s())) + 2) *
200           std::sqrt((std::exp(sqr(d.s())) - 1));
201     double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
202                       3*std::exp(2*sqr(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
test5()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(v);
223     }
224     double mean = std::accumulate(u.begin(), u.end(), 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(dbl);
232         var += d2;
233         skew += dbl * d2;
234         kurtosis += d2 * d2;
235     }
236     var /= u.size();
237     double dev = std::sqrt(var);
238     skew /= u.size() * dev * var;
239     kurtosis /= u.size() * var * var;
240     kurtosis -= 3;
241     double x_mean = std::exp(d.m() + sqr(d.s())/2);
242     double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
243     double x_skew = (std::exp(sqr(d.s())) + 2) *
244           std::sqrt((std::exp(sqr(d.s())) - 1));
245     double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
246                       3*std::exp(2*sqr(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 
main(int,char **)253 int main(int, char**)
254 {
255     test1();
256     test2();
257     test3();
258     test4();
259     test5();
260 
261   return 0;
262 }
263