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