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); |
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 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(x: 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 std::minstd_rand G; |
73 | G g; |
74 | D d(10); |
75 | const int N = 1000000; |
76 | std::vector<D::result_type> u; |
77 | for (int i = 0; i < N; ++i) |
78 | u.push_back(d(g)); |
79 | double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size(); |
80 | double var = 0; |
81 | double skew = 0; |
82 | double kurtosis = 0; |
83 | for (unsigned i = 0; i < u.size(); ++i) |
84 | { |
85 | double dbl = (u[i] - mean); |
86 | double d2 = sqr(dbl); |
87 | var += d2; |
88 | skew += dbl * d2; |
89 | kurtosis += d2 * d2; |
90 | } |
91 | var /= u.size(); |
92 | double dev = std::sqrt(x: var); |
93 | skew /= u.size() * dev * var; |
94 | kurtosis /= u.size() * var * var; |
95 | kurtosis -= 3; |
96 | double x_mean = 0; |
97 | double x_var = d.n() / (d.n() - 2); |
98 | double x_skew = 0; |
99 | double x_kurtosis = 6 / (d.n() - 4); |
100 | assert(std::abs(mean - x_mean) < 0.01); |
101 | assert(std::abs((var - x_var) / x_var) < 0.01); |
102 | assert(std::abs(skew - x_skew) < 0.01); |
103 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.04); |
104 | } |
105 | { |
106 | typedef std::student_t_distribution<> D; |
107 | typedef std::minstd_rand G; |
108 | G g; |
109 | D d(100); |
110 | const int N = 1000000; |
111 | std::vector<D::result_type> u; |
112 | for (int i = 0; i < N; ++i) |
113 | u.push_back(d(g)); |
114 | double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size(); |
115 | double var = 0; |
116 | double skew = 0; |
117 | double kurtosis = 0; |
118 | for (unsigned i = 0; i < u.size(); ++i) |
119 | { |
120 | double dbl = (u[i] - mean); |
121 | double d2 = sqr(dbl); |
122 | var += d2; |
123 | skew += dbl * d2; |
124 | kurtosis += d2 * d2; |
125 | } |
126 | var /= u.size(); |
127 | double dev = std::sqrt(x: var); |
128 | skew /= u.size() * dev * var; |
129 | kurtosis /= u.size() * var * var; |
130 | kurtosis -= 3; |
131 | double x_mean = 0; |
132 | double x_var = d.n() / (d.n() - 2); |
133 | double x_skew = 0; |
134 | double x_kurtosis = 6 / (d.n() - 4); |
135 | assert(std::abs(mean - x_mean) < 0.01); |
136 | assert(std::abs((var - x_var) / x_var) < 0.01); |
137 | assert(std::abs(skew - x_skew) < 0.01); |
138 | assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02); |
139 | } |
140 | |
141 | return 0; |
142 | } |
143 | |