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 weibull_distribution
16
17// template<class _URNG> result_type operator()(_URNG& g);
18
19#include <random>
20#include <cassert>
21#include <vector>
22#include <numeric>
23#include <cstddef>
24
25template <class T>
26inline
27T
28sqr(T x)
29{
30    return x * x;
31}
32
33int main()
34{
35    {
36        typedef std::weibull_distribution<> D;
37        typedef std::mt19937 G;
38        G g;
39        D d(0.5, 2);
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(d.min() <= v);
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 (std::size_t 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 = d.b() * std::tgamma(1 + 1/d.a());
66        double x_var = sqr(d.b()) * std::tgamma(1 + 2/d.a()) - sqr(x_mean);
67        double x_skew = (sqr(d.b())*d.b() * std::tgamma(1 + 3/d.a()) -
68                        3*x_mean*x_var - sqr(x_mean)*x_mean) /
69                        (std::sqrt(x_var)*x_var);
70        double x_kurtosis = (sqr(sqr(d.b())) * std::tgamma(1 + 4/d.a()) -
71                       4*x_skew*x_var*sqrt(x_var)*x_mean -
72                       6*sqr(x_mean)*x_var - sqr(sqr(x_mean))) / sqr(x_var) - 3;
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.01);
76        assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
77    }
78    {
79        typedef std::weibull_distribution<> D;
80        typedef std::mt19937 G;
81        G g;
82        D d(1, .5);
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);
88            assert(d.min() <= v);
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 (std::size_t i = 0; i < u.size(); ++i)
96        {
97            double dbl = (u[i] - mean);
98            double d2 = sqr(dbl);
99            var += d2;
100            skew += dbl * 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 = d.b() * std::tgamma(1 + 1/d.a());
109        double x_var = sqr(d.b()) * std::tgamma(1 + 2/d.a()) - sqr(x_mean);
110        double x_skew = (sqr(d.b())*d.b() * std::tgamma(1 + 3/d.a()) -
111                        3*x_mean*x_var - sqr(x_mean)*x_mean) /
112                        (std::sqrt(x_var)*x_var);
113        double x_kurtosis = (sqr(sqr(d.b())) * std::tgamma(1 + 4/d.a()) -
114                       4*x_skew*x_var*sqrt(x_var)*x_mean -
115                       6*sqr(x_mean)*x_var - sqr(sqr(x_mean))) / sqr(x_var) - 3;
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.01);
120    }
121    {
122        typedef std::weibull_distribution<> D;
123        typedef std::mt19937 G;
124        G g;
125        D d(2, 3);
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);
131            assert(d.min() <= v);
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 (std::size_t i = 0; i < u.size(); ++i)
139        {
140            double dbl = (u[i] - mean);
141            double d2 = sqr(dbl);
142            var += d2;
143            skew += dbl * 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 = d.b() * std::tgamma(1 + 1/d.a());
152        double x_var = sqr(d.b()) * std::tgamma(1 + 2/d.a()) - sqr(x_mean);
153        double x_skew = (sqr(d.b())*d.b() * std::tgamma(1 + 3/d.a()) -
154                        3*x_mean*x_var - sqr(x_mean)*x_mean) /
155                        (std::sqrt(x_var)*x_var);
156        double x_kurtosis = (sqr(sqr(d.b())) * std::tgamma(1 + 4/d.a()) -
157                       4*x_skew*x_var*sqrt(x_var)*x_mean -
158                       6*sqr(x_mean)*x_var - sqr(sqr(x_mean))) / sqr(x_var) - 3;
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.01);
162        assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
163    }
164}
165