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 student_t_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        typedef std::student_t_distribution<> D;
36        typedef D::param_type P;
37        typedef std::minstd_rand G;
38        G g;
39        D d;
40        P p(5.5);
41        const int N = 1000000;
42        std::vector<D::result_type> u;
43        for (int i = 0; i < N; ++i)
44            u.push_back(d(g, p));
45        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
46        double var = 0;
47        double skew = 0;
48        double kurtosis = 0;
49        for (unsigned i = 0; i < u.size(); ++i)
50        {
51            double dbl = (u[i] - mean);
52            double d2 = sqr(dbl);
53            var += d2;
54            skew += dbl * d2;
55            kurtosis += d2 * d2;
56        }
57        var /= u.size();
58        double dev = std::sqrt(var);
59        skew /= u.size() * dev * var;
60        kurtosis /= u.size() * var * var;
61        kurtosis -= 3;
62        double x_mean = 0;
63        double x_var = p.n() / (p.n() - 2);
64        double x_skew = 0;
65        double x_kurtosis = 6 / (p.n() - 4);
66        assert(std::abs(mean - x_mean) < 0.01);
67        assert(std::abs((var - x_var) / x_var) < 0.01);
68        assert(std::abs(skew - x_skew) < 0.01);
69        assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.2);
70    }
71    {
72        typedef std::student_t_distribution<> D;
73        typedef D::param_type P;
74        typedef std::minstd_rand G;
75        G g;
76        D d;
77        P p(10);
78        const int N = 1000000;
79        std::vector<D::result_type> u;
80        for (int i = 0; i < N; ++i)
81            u.push_back(d(g, p));
82        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
83        double var = 0;
84        double skew = 0;
85        double kurtosis = 0;
86        for (unsigned i = 0; i < u.size(); ++i)
87        {
88            double dbl = (u[i] - mean);
89            double d2 = sqr(dbl);
90            var += d2;
91            skew += dbl * d2;
92            kurtosis += d2 * d2;
93        }
94        var /= u.size();
95        double dev = std::sqrt(var);
96        skew /= u.size() * dev * var;
97        kurtosis /= u.size() * var * var;
98        kurtosis -= 3;
99        double x_mean = 0;
100        double x_var = p.n() / (p.n() - 2);
101        double x_skew = 0;
102        double x_kurtosis = 6 / (p.n() - 4);
103        assert(std::abs(mean - x_mean) < 0.01);
104        assert(std::abs((var - x_var) / x_var) < 0.01);
105        assert(std::abs(skew - x_skew) < 0.01);
106        assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.04);
107    }
108    {
109        typedef std::student_t_distribution<> D;
110        typedef D::param_type P;
111        typedef std::minstd_rand G;
112        G g;
113        D d;
114        P p(100);
115        const int N = 1000000;
116        std::vector<D::result_type> u;
117        for (int i = 0; i < N; ++i)
118            u.push_back(d(g, p));
119        double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
120        double var = 0;
121        double skew = 0;
122        double kurtosis = 0;
123        for (unsigned i = 0; i < u.size(); ++i)
124        {
125            double dbl = (u[i] - mean);
126            double d2 = sqr(dbl);
127            var += d2;
128            skew += dbl * d2;
129            kurtosis += d2 * d2;
130        }
131        var /= u.size();
132        double dev = std::sqrt(var);
133        skew /= u.size() * dev * var;
134        kurtosis /= u.size() * var * var;
135        kurtosis -= 3;
136        double x_mean = 0;
137        double x_var = p.n() / (p.n() - 2);
138        double x_skew = 0;
139        double x_kurtosis = 6 / (p.n() - 4);
140        assert(std::abs(mean - x_mean) < 0.01);
141        assert(std::abs((var - x_var) / x_var) < 0.01);
142        assert(std::abs(skew - x_skew) < 0.01);
143        assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
144    }
145}
146