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 student_t_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
sqr(T x)25 sqr(T x)
26 {
27 return x * x;
28 }
29
main()30 int main()
31 {
32 {
33 typedef std::student_t_distribution<> D;
34 typedef D::param_type P;
35 typedef std::minstd_rand G;
36 G g;
37 D d;
38 P p(5.5);
39 const int N = 1000000;
40 std::vector<D::result_type> u;
41 for (int i = 0; i < N; ++i)
42 u.push_back(d(g, p));
43 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
44 double var = 0;
45 double skew = 0;
46 double kurtosis = 0;
47 for (int i = 0; i < u.size(); ++i)
48 {
49 double d = (u[i] - mean);
50 double d2 = sqr(d);
51 var += d2;
52 skew += d * d2;
53 kurtosis += d2 * d2;
54 }
55 var /= u.size();
56 double dev = std::sqrt(var);
57 skew /= u.size() * dev * var;
58 kurtosis /= u.size() * var * var;
59 kurtosis -= 3;
60 double x_mean = 0;
61 double x_var = p.n() / (p.n() - 2);
62 double x_skew = 0;
63 double x_kurtosis = 6 / (p.n() - 4);
64 assert(std::abs(mean - x_mean) < 0.01);
65 assert(std::abs((var - x_var) / x_var) < 0.01);
66 assert(std::abs(skew - x_skew) < 0.01);
67 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.2);
68 }
69 {
70 typedef std::student_t_distribution<> D;
71 typedef D::param_type P;
72 typedef std::minstd_rand G;
73 G g;
74 D d;
75 P p(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, p));
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 (int i = 0; i < u.size(); ++i)
85 {
86 double d = (u[i] - mean);
87 double d2 = sqr(d);
88 var += d2;
89 skew += d * 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 = p.n() / (p.n() - 2);
99 double x_skew = 0;
100 double x_kurtosis = 6 / (p.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;
112 P p(100);
113 const int N = 1000000;
114 std::vector<D::result_type> u;
115 for (int i = 0; i < N; ++i)
116 u.push_back(d(g, p));
117 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
118 double var = 0;
119 double skew = 0;
120 double kurtosis = 0;
121 for (int i = 0; i < u.size(); ++i)
122 {
123 double d = (u[i] - mean);
124 double d2 = sqr(d);
125 var += d2;
126 skew += d * d2;
127 kurtosis += d2 * d2;
128 }
129 var /= u.size();
130 double dev = std::sqrt(var);
131 skew /= u.size() * dev * var;
132 kurtosis /= u.size() * var * var;
133 kurtosis -= 3;
134 double x_mean = 0;
135 double x_var = p.n() / (p.n() - 2);
136 double x_skew = 0;
137 double x_kurtosis = 6 / (p.n() - 4);
138 assert(std::abs(mean - x_mean) < 0.01);
139 assert(std::abs((var - x_var) / x_var) < 0.01);
140 assert(std::abs(skew - x_skew) < 0.01);
141 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
142 }
143 }
144