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