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 // <random> 11 12 // template<class IntType = int> 13 // class poisson_distribution 14 15 // template<class _URNG> result_type operator()(_URNG& g); 16 17 #include <random> 18 #include <cassert> 19 #include <vector> 20 #include <numeric> 21 22 template <class T> 23 inline 24 T sqr(T x)25sqr(T x) 26 { 27 return x * x; 28 } 29 main()30int main() 31 { 32 { 33 typedef std::poisson_distribution<> D; 34 typedef std::minstd_rand G; 35 G g; 36 D d(2); 37 const int N = 100000; 38 std::vector<double> u; 39 for (int i = 0; i < N; ++i) 40 { 41 D::result_type v = d(g); 42 assert(d.min() <= v && v <= d.max()); 43 u.push_back(v); 44 } 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 (int i = 0; i < u.size(); ++i) 50 { 51 double d = (u[i] - mean); 52 double d2 = sqr(d); 53 var += d2; 54 skew += d * 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 = d.mean(); 63 double x_var = d.mean(); 64 double x_skew = 1 / std::sqrt(x_var); 65 double x_kurtosis = 1 / x_var; 66 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 67 assert(std::abs((var - x_var) / x_var) < 0.01); 68 assert(std::abs((skew - x_skew) / x_skew) < 0.01); 69 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03); 70 } 71 { 72 typedef std::poisson_distribution<> D; 73 typedef std::minstd_rand G; 74 G g; 75 D d(0.75); 76 const int N = 100000; 77 std::vector<double> u; 78 for (int i = 0; i < N; ++i) 79 { 80 D::result_type v = d(g); 81 assert(d.min() <= v && v <= d.max()); 82 u.push_back(v); 83 } 84 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size(); 85 double var = 0; 86 double skew = 0; 87 double kurtosis = 0; 88 for (int i = 0; i < u.size(); ++i) 89 { 90 double d = (u[i] - mean); 91 double d2 = sqr(d); 92 var += d2; 93 skew += d * d2; 94 kurtosis += d2 * d2; 95 } 96 var /= u.size(); 97 double dev = std::sqrt(var); 98 skew /= u.size() * dev * var; 99 kurtosis /= u.size() * var * var; 100 kurtosis -= 3; 101 double x_mean = d.mean(); 102 double x_var = d.mean(); 103 double x_skew = 1 / std::sqrt(x_var); 104 double x_kurtosis = 1 / x_var; 105 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 106 assert(std::abs((var - x_var) / x_var) < 0.01); 107 assert(std::abs((skew - x_skew) / x_skew) < 0.01); 108 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.04); 109 } 110 { 111 typedef std::poisson_distribution<> D; 112 typedef std::mt19937 G; 113 G g; 114 D d(20); 115 const int N = 1000000; 116 std::vector<double> u; 117 for (int i = 0; i < N; ++i) 118 { 119 D::result_type v = d(g); 120 assert(d.min() <= v && v <= d.max()); 121 u.push_back(v); 122 } 123 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size(); 124 double var = 0; 125 double skew = 0; 126 double kurtosis = 0; 127 for (int i = 0; i < u.size(); ++i) 128 { 129 double d = (u[i] - mean); 130 double d2 = sqr(d); 131 var += d2; 132 skew += d * d2; 133 kurtosis += d2 * d2; 134 } 135 var /= u.size(); 136 double dev = std::sqrt(var); 137 skew /= u.size() * dev * var; 138 kurtosis /= u.size() * var * var; 139 kurtosis -= 3; 140 double x_mean = d.mean(); 141 double x_var = d.mean(); 142 double x_skew = 1 / std::sqrt(x_var); 143 double x_kurtosis = 1 / x_var; 144 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 145 assert(std::abs((var - x_var) / x_var) < 0.01); 146 assert(std::abs((skew - x_skew) / x_skew) < 0.01); 147 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 148 } 149 } 150