1f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant//===----------------------------------------------------------------------===// 2f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant// 3f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant// The LLVM Compiler Infrastructure 4f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant// 5b64f8b07c104c6cc986570ac8ee0ed16a9f23976Howard Hinnant// This file is dual licensed under the MIT and the University of Illinois Open 6b64f8b07c104c6cc986570ac8ee0ed16a9f23976Howard Hinnant// Source Licenses. See LICENSE.TXT for details. 7f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant// 8f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant//===----------------------------------------------------------------------===// 9f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant 10f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant// <random> 11f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant 12f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant// template<class IntType = int> 13f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant// class negative_binomial_distribution 14f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant 15f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant// template<class _URNG> result_type operator()(_URNG& g, const param_type& parm); 16f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant 17f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant#include <random> 18f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant#include <numeric> 19f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant#include <vector> 20f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant#include <cassert> 21f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant 22f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnanttemplate <class T> 23f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnantinline 24f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward HinnantT 25f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnantsqr(T x) 26f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant{ 27f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant return x * x; 28f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant} 29f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant 30f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnantint main() 31f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant{ 32f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant { 33f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant typedef std::negative_binomial_distribution<> D; 34f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant typedef D::param_type P; 35f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant typedef std::minstd_rand G; 36f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant G g; 37f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant D d(16, .75); 38f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant P p(5, .75); 39f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant const int N = 1000000; 40f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant std::vector<D::result_type> u; 41f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant for (int i = 0; i < N; ++i) 42f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant { 43f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant D::result_type v = d(g, p); 44f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant assert(d.min() <= v && v <= d.max()); 45f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant u.push_back(v); 46f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant } 47f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double mean = std::accumulate(u.begin(), u.end(), 48f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double(0)) / u.size(); 49f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double var = 0; 50f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double skew = 0; 51f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double kurtosis = 0; 52f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant for (int i = 0; i < u.size(); ++i) 53f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant { 54f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double d = (u[i] - mean); 55f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double d2 = sqr(d); 56f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant var += d2; 57f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant skew += d * d2; 58f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant kurtosis += d2 * d2; 59f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant } 60f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant var /= u.size(); 61f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double dev = std::sqrt(var); 62f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant skew /= u.size() * dev * var; 63f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant kurtosis /= u.size() * var * var; 64f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant kurtosis -= 3; 65f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double x_mean = p.k() * (1 - p.p()) / p.p(); 66f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double x_var = x_mean / p.p(); 67f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double x_skew = (2 - p.p()) / std::sqrt(p.k() * (1 - p.p())); 68f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double x_kurtosis = 6. / p.k() + sqr(p.p()) / (p.k() * (1 - p.p())); 69d6d1171f2c3f254582ae1d5b9e14cea0ea8e701bHoward Hinnant assert(std::abs((mean - x_mean) / x_mean) < 0.01); 70d6d1171f2c3f254582ae1d5b9e14cea0ea8e701bHoward Hinnant assert(std::abs((var - x_var) / x_var) < 0.01); 71d6d1171f2c3f254582ae1d5b9e14cea0ea8e701bHoward Hinnant assert(std::abs((skew - x_skew) / x_skew) < 0.01); 72d6d1171f2c3f254582ae1d5b9e14cea0ea8e701bHoward Hinnant assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 73f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant } 74f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant { 75f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant typedef std::negative_binomial_distribution<> D; 76f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant typedef D::param_type P; 77f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant typedef std::mt19937 G; 78f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant G g; 79f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant D d(16, .75); 80f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant P p(30, .03125); 81f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant const int N = 1000000; 82f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant std::vector<D::result_type> u; 83f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant for (int i = 0; i < N; ++i) 84f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant { 85f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant D::result_type v = d(g, p); 86f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant assert(d.min() <= v && v <= d.max()); 87f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant u.push_back(v); 88f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant } 89f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double mean = std::accumulate(u.begin(), u.end(), 90f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double(0)) / u.size(); 91f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double var = 0; 92f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double skew = 0; 93f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double kurtosis = 0; 94f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant for (int i = 0; i < u.size(); ++i) 95f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant { 96f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double d = (u[i] - mean); 97f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double d2 = sqr(d); 98f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant var += d2; 99f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant skew += d * d2; 100f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant kurtosis += d2 * d2; 101f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant } 102f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant var /= u.size(); 103f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double dev = std::sqrt(var); 104f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant skew /= u.size() * dev * var; 105f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant kurtosis /= u.size() * var * var; 106f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant kurtosis -= 3; 107f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double x_mean = p.k() * (1 - p.p()) / p.p(); 108f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double x_var = x_mean / p.p(); 109f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double x_skew = (2 - p.p()) / std::sqrt(p.k() * (1 - p.p())); 110f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double x_kurtosis = 6. / p.k() + sqr(p.p()) / (p.k() * (1 - p.p())); 111d6d1171f2c3f254582ae1d5b9e14cea0ea8e701bHoward Hinnant assert(std::abs((mean - x_mean) / x_mean) < 0.01); 112d6d1171f2c3f254582ae1d5b9e14cea0ea8e701bHoward Hinnant assert(std::abs((var - x_var) / x_var) < 0.01); 113d6d1171f2c3f254582ae1d5b9e14cea0ea8e701bHoward Hinnant assert(std::abs((skew - x_skew) / x_skew) < 0.01); 114d6d1171f2c3f254582ae1d5b9e14cea0ea8e701bHoward Hinnant assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01); 115f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant } 116f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant { 117f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant typedef std::negative_binomial_distribution<> D; 118f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant typedef D::param_type P; 119f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant typedef std::mt19937 G; 120f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant G g; 121f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant D d(16, .75); 122f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant P p(40, .25); 123f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant const int N = 1000000; 124f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant std::vector<D::result_type> u; 125f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant for (int i = 0; i < N; ++i) 126f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant { 127f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant D::result_type v = d(g, p); 128f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant assert(d.min() <= v && v <= d.max()); 129f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant u.push_back(v); 130f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant } 131f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double mean = std::accumulate(u.begin(), u.end(), 132f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double(0)) / u.size(); 133f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double var = 0; 134f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double skew = 0; 135f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double kurtosis = 0; 136f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant for (int i = 0; i < u.size(); ++i) 137f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant { 138f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double d = (u[i] - mean); 139f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double d2 = sqr(d); 140f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant var += d2; 141f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant skew += d * d2; 142f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant kurtosis += d2 * d2; 143f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant } 144f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant var /= u.size(); 145f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double dev = std::sqrt(var); 146f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant skew /= u.size() * dev * var; 147f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant kurtosis /= u.size() * var * var; 148f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant kurtosis -= 3; 149f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double x_mean = p.k() * (1 - p.p()) / p.p(); 150f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double x_var = x_mean / p.p(); 151f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double x_skew = (2 - p.p()) / std::sqrt(p.k() * (1 - p.p())); 152f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant double x_kurtosis = 6. / p.k() + sqr(p.p()) / (p.k() * (1 - p.p())); 153d6d1171f2c3f254582ae1d5b9e14cea0ea8e701bHoward Hinnant assert(std::abs((mean - x_mean) / x_mean) < 0.01); 154d6d1171f2c3f254582ae1d5b9e14cea0ea8e701bHoward Hinnant assert(std::abs((var - x_var) / x_var) < 0.01); 155d6d1171f2c3f254582ae1d5b9e14cea0ea8e701bHoward Hinnant assert(std::abs((skew - x_skew) / x_skew) < 0.01); 156d6d1171f2c3f254582ae1d5b9e14cea0ea8e701bHoward Hinnant assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03); 157f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant } 158f2fe5d5a61d1f5842068793620aaff1e88173f7fHoward Hinnant} 159