1 //===----------------------------------------------------------------------===//
2 //
3 //                     The LLVM Compiler Infrastructure
4 //
5 // This file is dual licensed under the MIT and the University of Illinois Open
6 // Source Licenses. See LICENSE.TXT for details.
7 //
8 //===----------------------------------------------------------------------===//
9 
10 // <random>
11 
12 // template<class RealType = double>
13 // class lognormal_distribution
14 
15 // template<class _URNG> result_type operator()(_URNG& g);
16 
17 #include <random>
18 #include <cassert>
19 #include <vector>
20 #include <numeric>
21 
22 template <class T>
23 inline
24 T
25 sqr(T x)
26 {
27     return x * x;
28 }
29 
30 int main()
31 {
32     {
33         typedef std::lognormal_distribution<> D;
34         typedef D::param_type P;
35         typedef std::mt19937 G;
36         G g;
37         D d(-1./8192, 0.015625);
38         const int N = 1000000;
39         std::vector<D::result_type> u;
40         for (int i = 0; i < N; ++i)
41         {
42             D::result_type v = d(g);
43             assert(v > 0);
44             u.push_back(v);
45         }
46         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
47         double var = 0;
48         double skew = 0;
49         double kurtosis = 0;
50         for (int i = 0; i < u.size(); ++i)
51         {
52             double d = (u[i] - mean);
53             double d2 = sqr(d);
54             var += d2;
55             skew += d * d2;
56             kurtosis += d2 * d2;
57         }
58         var /= u.size();
59         double dev = std::sqrt(var);
60         skew /= u.size() * dev * var;
61         kurtosis /= u.size() * var * var;
62         kurtosis -= 3;
63         double x_mean = std::exp(d.m() + sqr(d.s())/2);
64         double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
65         double x_skew = (std::exp(sqr(d.s())) + 2) *
66               std::sqrt((std::exp(sqr(d.s())) - 1));
67         double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
68                           3*std::exp(2*sqr(d.s())) - 6;
69         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
70         assert(std::abs((var - x_var) / x_var) < 0.01);
71         assert(std::abs((skew - x_skew) / x_skew) < 0.05);
72         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.25);
73     }
74     {
75         typedef std::lognormal_distribution<> D;
76         typedef D::param_type P;
77         typedef std::mt19937 G;
78         G g;
79         D d(-1./32, 0.25);
80         const int N = 1000000;
81         std::vector<D::result_type> u;
82         for (int i = 0; i < N; ++i)
83         {
84             D::result_type v = d(g);
85             assert(v > 0);
86             u.push_back(v);
87         }
88         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
89         double var = 0;
90         double skew = 0;
91         double kurtosis = 0;
92         for (int i = 0; i < u.size(); ++i)
93         {
94             double d = (u[i] - mean);
95             double d2 = sqr(d);
96             var += d2;
97             skew += d * d2;
98             kurtosis += d2 * d2;
99         }
100         var /= u.size();
101         double dev = std::sqrt(var);
102         skew /= u.size() * dev * var;
103         kurtosis /= u.size() * var * var;
104         kurtosis -= 3;
105         double x_mean = std::exp(d.m() + sqr(d.s())/2);
106         double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
107         double x_skew = (std::exp(sqr(d.s())) + 2) *
108               std::sqrt((std::exp(sqr(d.s())) - 1));
109         double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
110                           3*std::exp(2*sqr(d.s())) - 6;
111         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
112         assert(std::abs((var - x_var) / x_var) < 0.01);
113         assert(std::abs((skew - x_skew) / x_skew) < 0.01);
114         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
115     }
116     {
117         typedef std::lognormal_distribution<> D;
118         typedef D::param_type P;
119         typedef std::mt19937 G;
120         G g;
121         D d(-1./8, 0.5);
122         const int N = 1000000;
123         std::vector<D::result_type> u;
124         for (int i = 0; i < N; ++i)
125         {
126             D::result_type v = d(g);
127             assert(v > 0);
128             u.push_back(v);
129         }
130         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
131         double var = 0;
132         double skew = 0;
133         double kurtosis = 0;
134         for (int i = 0; i < u.size(); ++i)
135         {
136             double d = (u[i] - mean);
137             double d2 = sqr(d);
138             var += d2;
139             skew += d * d2;
140             kurtosis += d2 * d2;
141         }
142         var /= u.size();
143         double dev = std::sqrt(var);
144         skew /= u.size() * dev * var;
145         kurtosis /= u.size() * var * var;
146         kurtosis -= 3;
147         double x_mean = std::exp(d.m() + sqr(d.s())/2);
148         double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
149         double x_skew = (std::exp(sqr(d.s())) + 2) *
150               std::sqrt((std::exp(sqr(d.s())) - 1));
151         double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
152                           3*std::exp(2*sqr(d.s())) - 6;
153         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
154         assert(std::abs((var - x_var) / x_var) < 0.01);
155         assert(std::abs((skew - x_skew) / x_skew) < 0.02);
156         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.05);
157     }
158     {
159         typedef std::lognormal_distribution<> D;
160         typedef D::param_type P;
161         typedef std::mt19937 G;
162         G g;
163         D d;
164         const int N = 1000000;
165         std::vector<D::result_type> u;
166         for (int i = 0; i < N; ++i)
167         {
168             D::result_type v = d(g);
169             assert(v > 0);
170             u.push_back(v);
171         }
172         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
173         double var = 0;
174         double skew = 0;
175         double kurtosis = 0;
176         for (int i = 0; i < u.size(); ++i)
177         {
178             double d = (u[i] - mean);
179             double d2 = sqr(d);
180             var += d2;
181             skew += d * d2;
182             kurtosis += d2 * d2;
183         }
184         var /= u.size();
185         double dev = std::sqrt(var);
186         skew /= u.size() * dev * var;
187         kurtosis /= u.size() * var * var;
188         kurtosis -= 3;
189         double x_mean = std::exp(d.m() + sqr(d.s())/2);
190         double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
191         double x_skew = (std::exp(sqr(d.s())) + 2) *
192               std::sqrt((std::exp(sqr(d.s())) - 1));
193         double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
194                           3*std::exp(2*sqr(d.s())) - 6;
195         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
196         assert(std::abs((var - x_var) / x_var) < 0.02);
197         assert(std::abs((skew - x_skew) / x_skew) < 0.08);
198         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.4);
199     }
200     {
201         typedef std::lognormal_distribution<> D;
202         typedef D::param_type P;
203         typedef std::mt19937 G;
204         G g;
205         D d(-0.78125, 1.25);
206         const int N = 1000000;
207         std::vector<D::result_type> u;
208         for (int i = 0; i < N; ++i)
209         {
210             D::result_type v = d(g);
211             assert(v > 0);
212             u.push_back(v);
213         }
214         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
215         double var = 0;
216         double skew = 0;
217         double kurtosis = 0;
218         for (int i = 0; i < u.size(); ++i)
219         {
220             double d = (u[i] - mean);
221             double d2 = sqr(d);
222             var += d2;
223             skew += d * d2;
224             kurtosis += d2 * d2;
225         }
226         var /= u.size();
227         double dev = std::sqrt(var);
228         skew /= u.size() * dev * var;
229         kurtosis /= u.size() * var * var;
230         kurtosis -= 3;
231         double x_mean = std::exp(d.m() + sqr(d.s())/2);
232         double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
233         double x_skew = (std::exp(sqr(d.s())) + 2) *
234               std::sqrt((std::exp(sqr(d.s())) - 1));
235         double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
236                           3*std::exp(2*sqr(d.s())) - 6;
237         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
238         assert(std::abs((var - x_var) / x_var) < 0.04);
239         assert(std::abs((skew - x_skew) / x_skew) < 0.2);
240         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.7);
241     }
242 }
243