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