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// REQUIRES: long_tests
11
12// <random>
13
14// template<class RealType = double>
15// class lognormal_distribution
16
17// template<class _URNG> result_type operator()(_URNG& g, const param_type& parm);
18
19#include <random>
20#include <cassert>
21#include <vector>
22#include <numeric>
23
24template <class T>
25inline
26T
27sqr(T x)
28{
29    return x * x;
30}
31
32int main()
33{
34
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 (int i = 0; i < u.size(); ++i)
55        {
56            double d = (u[i] - mean);
57            double d2 = sqr(d);
58            var += d2;
59            skew += d * 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        typedef std::lognormal_distribution<> D;
80        typedef D::param_type P;
81        typedef std::mt19937 G;
82        G g;
83        D d;
84        P p(-1./32, 0.25);
85        const int N = 1000000;
86        std::vector<D::result_type> u;
87        for (int i = 0; i < N; ++i)
88        {
89            D::result_type v = d(g, p);
90            assert(v > 0);
91            u.push_back(v);
92        }
93        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
94        double var = 0;
95        double skew = 0;
96        double kurtosis = 0;
97        for (int i = 0; i < u.size(); ++i)
98        {
99            double d = (u[i] - mean);
100            double d2 = sqr(d);
101            var += d2;
102            skew += d * d2;
103            kurtosis += d2 * d2;
104        }
105        var /= u.size();
106        double dev = std::sqrt(var);
107        skew /= u.size() * dev * var;
108        kurtosis /= u.size() * var * var;
109        kurtosis -= 3;
110        double x_mean = std::exp(p.m() + sqr(p.s())/2);
111        double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
112        double x_skew = (std::exp(sqr(p.s())) + 2) *
113              std::sqrt((std::exp(sqr(p.s())) - 1));
114        double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
115                          3*std::exp(2*sqr(p.s())) - 6;
116        assert(std::abs((mean - x_mean) / x_mean) < 0.01);
117        assert(std::abs((var - x_var) / x_var) < 0.01);
118        assert(std::abs((skew - x_skew) / x_skew) < 0.01);
119        assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
120    }
121    {
122        typedef std::lognormal_distribution<> D;
123        typedef D::param_type P;
124        typedef std::mt19937 G;
125        G g;
126        D d;
127        P p(-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, p);
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 (int i = 0; i < u.size(); ++i)
141        {
142            double d = (u[i] - mean);
143            double d2 = sqr(d);
144            var += d2;
145            skew += d * 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(p.m() + sqr(p.s())/2);
154        double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
155        double x_skew = (std::exp(sqr(p.s())) + 2) *
156              std::sqrt((std::exp(sqr(p.s())) - 1));
157        double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
158                          3*std::exp(2*sqr(p.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        typedef std::lognormal_distribution<> D;
166        typedef D::param_type P;
167        typedef std::mt19937 G;
168        G g;
169        D d(3, 4);
170        P p;
171        const int N = 1000000;
172        std::vector<D::result_type> u;
173        for (int i = 0; i < N; ++i)
174        {
175            D::result_type v = d(g, p);
176            assert(v > 0);
177            u.push_back(v);
178        }
179        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
180        double var = 0;
181        double skew = 0;
182        double kurtosis = 0;
183        for (int i = 0; i < u.size(); ++i)
184        {
185            double d = (u[i] - mean);
186            double d2 = sqr(d);
187            var += d2;
188            skew += d * d2;
189            kurtosis += d2 * d2;
190        }
191        var /= u.size();
192        double dev = std::sqrt(var);
193        skew /= u.size() * dev * var;
194        kurtosis /= u.size() * var * var;
195        kurtosis -= 3;
196        double x_mean = std::exp(p.m() + sqr(p.s())/2);
197        double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
198        double x_skew = (std::exp(sqr(p.s())) + 2) *
199              std::sqrt((std::exp(sqr(p.s())) - 1));
200        double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
201                          3*std::exp(2*sqr(p.s())) - 6;
202        assert(std::abs((mean - x_mean) / x_mean) < 0.01);
203        assert(std::abs((var - x_var) / x_var) < 0.02);
204        assert(std::abs((skew - x_skew) / x_skew) < 0.08);
205        assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.4);
206    }
207    {
208        typedef std::lognormal_distribution<> D;
209        typedef D::param_type P;
210        typedef std::mt19937 G;
211        G g;
212        D d;
213        P p(-0.78125, 1.25);
214        const int N = 1000000;
215        std::vector<D::result_type> u;
216        for (int i = 0; i < N; ++i)
217        {
218            D::result_type v = d(g, p);
219            assert(v > 0);
220            u.push_back(v);
221        }
222        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
223        double var = 0;
224        double skew = 0;
225        double kurtosis = 0;
226        for (int i = 0; i < u.size(); ++i)
227        {
228            double d = (u[i] - mean);
229            double d2 = sqr(d);
230            var += d2;
231            skew += d * d2;
232            kurtosis += d2 * d2;
233        }
234        var /= u.size();
235        double dev = std::sqrt(var);
236        skew /= u.size() * dev * var;
237        kurtosis /= u.size() * var * var;
238        kurtosis -= 3;
239        double x_mean = std::exp(p.m() + sqr(p.s())/2);
240        double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
241        double x_skew = (std::exp(sqr(p.s())) + 2) *
242              std::sqrt((std::exp(sqr(p.s())) - 1));
243        double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
244                          3*std::exp(2*sqr(p.s())) - 6;
245        assert(std::abs((mean - x_mean) / x_mean) < 0.01);
246        assert(std::abs((var - x_var) / x_var) < 0.04);
247        assert(std::abs((skew - x_skew) / x_skew) < 0.2);
248        assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.7);
249    }
250}
251