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 IntType = int> 14 // class binomial_distribution 15 16 // template<class _URNG> result_type operator()(_URNG& g, const param_type& parm); 17 18 #include <random> 19 #include <numeric> 20 #include <vector> 21 #include <cassert> 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::binomial_distribution<> D; 37 typedef D::param_type P; 38 typedef std::mt19937_64 G; 39 G g; 40 D d(16, .75); 41 P p(5, .75); 42 const int N = 1000000; 43 std::vector<D::result_type> u; 44 for (int i = 0; i < N; ++i) 45 { 46 D::result_type v = d(g, p); 47 assert(0 <= v && v <= p.t()); 48 u.push_back(v); 49 } 50 double mean = std::accumulate(u.begin(), u.end(), 51 double(0)) / u.size(); 52 double var = 0; 53 double skew = 0; 54 double kurtosis = 0; 55 for (unsigned i = 0; i < u.size(); ++i) 56 { 57 double dbl = (u[i] - mean); 58 double d2 = sqr(dbl); 59 var += d2; 60 skew += dbl * d2; 61 kurtosis += d2 * d2; 62 } 63 var /= u.size(); 64 double dev = std::sqrt(var); 65 skew /= u.size() * dev * var; 66 kurtosis /= u.size() * var * var; 67 kurtosis -= 3; 68 double x_mean = p.t() * p.p(); 69 double x_var = x_mean*(1-p.p()); 70 double x_skew = (1-2*p.p()) / std::sqrt(x_var); 71 double x_kurtosis = (1-6*p.p()*(1-p.p())) / x_var; 72 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 73 assert(std::abs((var - x_var) / x_var) < 0.01); 74 assert(std::abs((skew - x_skew) / x_skew) < 0.01); 75 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.04); 76 } 77 { 78 typedef std::binomial_distribution<> D; 79 typedef D::param_type P; 80 typedef std::mt19937 G; 81 G g; 82 D d(16, .75); 83 P p(30, .03125); 84 const int N = 100000; 85 std::vector<D::result_type> u; 86 for (int i = 0; i < N; ++i) 87 { 88 D::result_type v = d(g, p); 89 assert(0 <= v && v <= p.t()); 90 u.push_back(v); 91 } 92 double mean = std::accumulate(u.begin(), u.end(), 93 double(0)) / u.size(); 94 double var = 0; 95 double skew = 0; 96 double kurtosis = 0; 97 for (unsigned i = 0; i < u.size(); ++i) 98 { 99 double dbl = (u[i] - mean); 100 double d2 = sqr(dbl); 101 var += d2; 102 skew += dbl * d2; 103 kurtosis += d2 * d2; 104 } 105 var /= u.size(); 106 double dev = std::sqrt(var); 107 skew /= u.size() * dev * var; 108 kurtosis /= u.size() * var * var; 109 kurtosis -= 3; 110 double x_mean = p.t() * p.p(); 111 double x_var = x_mean*(1-p.p()); 112 double x_skew = (1-2*p.p()) / std::sqrt(x_var); 113 double x_kurtosis = (1-6*p.p()*(1-p.p())) / x_var; 114 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 115 assert(std::abs((var - x_var) / x_var) < 0.01); 116 assert(std::abs((skew - x_skew) / x_skew) < 0.01); 117 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 118 } 119 { 120 typedef std::binomial_distribution<> D; 121 typedef D::param_type P; 122 typedef std::mt19937 G; 123 G g; 124 D d(16, .75); 125 P p(40, .25); 126 const int N = 1000000; 127 std::vector<D::result_type> u; 128 for (int i = 0; i < N; ++i) 129 { 130 D::result_type v = d(g, p); 131 assert(0 <= v && v <= p.t()); 132 u.push_back(v); 133 } 134 double mean = std::accumulate(u.begin(), u.end(), 135 double(0)) / u.size(); 136 double var = 0; 137 double skew = 0; 138 double kurtosis = 0; 139 for (unsigned i = 0; i < u.size(); ++i) 140 { 141 double dbl = (u[i] - mean); 142 double d2 = sqr(dbl); 143 var += d2; 144 skew += dbl * d2; 145 kurtosis += d2 * d2; 146 } 147 var /= u.size(); 148 double dev = std::sqrt(var); 149 skew /= u.size() * dev * var; 150 kurtosis /= u.size() * var * var; 151 kurtosis -= 3; 152 double x_mean = p.t() * p.p(); 153 double x_var = x_mean*(1-p.p()); 154 double x_skew = (1-2*p.p()) / std::sqrt(x_var); 155 double x_kurtosis = (1-6*p.p()*(1-p.p())) / x_var; 156 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 157 assert(std::abs((var - x_var) / x_var) < 0.01); 158 assert(std::abs((skew - x_skew) / x_skew) < 0.04); 159 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.3); 160 } 161 162 return 0; 163 } 164