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