1 | //===----------------------------------------------------------------------===// |
2 | // |
3 | // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. |
4 | // See https://llvm.org/LICENSE.txt for license information. |
5 | // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception |
6 | // |
7 | //===----------------------------------------------------------------------===// |
8 | // |
9 | // REQUIRES: long_tests |
10 | |
11 | // <random> |
12 | |
13 | // template<class RealType = double> |
14 | // class student_t_distribution |
15 | |
16 | // template<class _URNG> result_type operator()(_URNG& g, const param_type& parm); |
17 | |
18 | #include <random> |
19 | #include <cassert> |
20 | #include <vector> |
21 | #include <numeric> |
22 | |
23 | #include "test_macros.h" |
24 | |
25 | template <class T> |
26 | inline |
27 | T |
28 | sqr(T x) |
29 | { |
30 | return x * x; |
31 | } |
32 | |
33 | int main(int, char**) |
34 | { |
35 | { |
36 | typedef std::student_t_distribution<> D; |
37 | typedef D::param_type P; |
38 | typedef std::minstd_rand G; |
39 | G g; |
40 | D d; |
41 | P p(5.5); |
42 | const int N = 1000000; |
43 | std::vector<D::result_type> u; |
44 | for (int i = 0; i < N; ++i) |
45 | u.push_back(d(g, p)); |
46 | double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size(); |
47 | double var = 0; |
48 | double skew = 0; |
49 | double kurtosis = 0; |
50 | for (unsigned i = 0; i < u.size(); ++i) |
51 | { |
52 | double dbl = (u[i] - mean); |
53 | double d2 = sqr(dbl); |
54 | var += d2; |
55 | skew += dbl * d2; |
56 | kurtosis += d2 * d2; |
57 | } |
58 | var /= u.size(); |
59 | double dev = std::sqrt(x: var); |
60 | skew /= u.size() * dev * var; |
61 | kurtosis /= u.size() * var * var; |
62 | kurtosis -= 3; |
63 | double x_mean = 0; |
64 | double x_var = p.n() / (p.n() - 2); |
65 | double x_skew = 0; |
66 | double x_kurtosis = 6 / (p.n() - 4); |
67 | assert(std::abs(mean - x_mean) < 0.01); |
68 | assert(std::abs((var - x_var) / x_var) < 0.01); |
69 | assert(std::abs(skew - x_skew) < 0.01); |
70 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.2); |
71 | } |
72 | { |
73 | typedef std::student_t_distribution<> D; |
74 | typedef D::param_type P; |
75 | typedef std::minstd_rand G; |
76 | G g; |
77 | D d; |
78 | P p(10); |
79 | const int N = 1000000; |
80 | std::vector<D::result_type> u; |
81 | for (int i = 0; i < N; ++i) |
82 | u.push_back(d(g, p)); |
83 | double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size(); |
84 | double var = 0; |
85 | double skew = 0; |
86 | double kurtosis = 0; |
87 | for (unsigned i = 0; i < u.size(); ++i) |
88 | { |
89 | double dbl = (u[i] - mean); |
90 | double d2 = sqr(dbl); |
91 | var += d2; |
92 | skew += dbl * d2; |
93 | kurtosis += d2 * d2; |
94 | } |
95 | var /= u.size(); |
96 | double dev = std::sqrt(x: var); |
97 | skew /= u.size() * dev * var; |
98 | kurtosis /= u.size() * var * var; |
99 | kurtosis -= 3; |
100 | double x_mean = 0; |
101 | double x_var = p.n() / (p.n() - 2); |
102 | double x_skew = 0; |
103 | double x_kurtosis = 6 / (p.n() - 4); |
104 | assert(std::abs(mean - x_mean) < 0.01); |
105 | assert(std::abs((var - x_var) / x_var) < 0.01); |
106 | assert(std::abs(skew - x_skew) < 0.01); |
107 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.04); |
108 | } |
109 | { |
110 | typedef std::student_t_distribution<> D; |
111 | typedef D::param_type P; |
112 | typedef std::minstd_rand G; |
113 | G g; |
114 | D d; |
115 | P p(100); |
116 | const int N = 1000000; |
117 | std::vector<D::result_type> u; |
118 | for (int i = 0; i < N; ++i) |
119 | u.push_back(d(g, p)); |
120 | double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size(); |
121 | double var = 0; |
122 | double skew = 0; |
123 | double kurtosis = 0; |
124 | for (unsigned i = 0; i < u.size(); ++i) |
125 | { |
126 | double dbl = (u[i] - mean); |
127 | double d2 = sqr(dbl); |
128 | var += d2; |
129 | skew += dbl * d2; |
130 | kurtosis += d2 * d2; |
131 | } |
132 | var /= u.size(); |
133 | double dev = std::sqrt(x: var); |
134 | skew /= u.size() * dev * var; |
135 | kurtosis /= u.size() * var * var; |
136 | kurtosis -= 3; |
137 | double x_mean = 0; |
138 | double x_var = p.n() / (p.n() - 2); |
139 | double x_skew = 0; |
140 | double x_kurtosis = 6 / (p.n() - 4); |
141 | assert(std::abs(mean - x_mean) < 0.01); |
142 | assert(std::abs((var - x_var) / x_var) < 0.01); |
143 | assert(std::abs(skew - x_skew) < 0.01); |
144 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02); |
145 | } |
146 | |
147 | return 0; |
148 | } |
149 | |