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