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// REQUIRES: long_tests 11 12// <random> 13 14// template<class RealType = double> 15// class lognormal_distribution 16 17// template<class _URNG> result_type operator()(_URNG& g, const param_type& parm); 18 19#include <random> 20#include <cassert> 21#include <vector> 22#include <numeric> 23 24template <class T> 25inline 26T 27sqr(T x) 28{ 29 return x * x; 30} 31 32int main() 33{ 34 35 { 36 typedef std::lognormal_distribution<> D; 37 typedef D::param_type P; 38 typedef std::mt19937 G; 39 G g; 40 D d; 41 P p(-1./8192, 0.015625); 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(v > 0); 48 u.push_back(v); 49 } 50 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size(); 51 double var = 0; 52 double skew = 0; 53 double kurtosis = 0; 54 for (int i = 0; i < u.size(); ++i) 55 { 56 double d = (u[i] - mean); 57 double d2 = sqr(d); 58 var += d2; 59 skew += d * d2; 60 kurtosis += d2 * d2; 61 } 62 var /= u.size(); 63 double dev = std::sqrt(var); 64 skew /= u.size() * dev * var; 65 kurtosis /= u.size() * var * var; 66 kurtosis -= 3; 67 double x_mean = std::exp(p.m() + sqr(p.s())/2); 68 double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s())); 69 double x_skew = (std::exp(sqr(p.s())) + 2) * 70 std::sqrt((std::exp(sqr(p.s())) - 1)); 71 double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) + 72 3*std::exp(2*sqr(p.s())) - 6; 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.05); 76 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.25); 77 } 78 { 79 typedef std::lognormal_distribution<> D; 80 typedef D::param_type P; 81 typedef std::mt19937 G; 82 G g; 83 D d; 84 P p(-1./32, 0.25); 85 const int N = 1000000; 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(v > 0); 91 u.push_back(v); 92 } 93 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size(); 94 double var = 0; 95 double skew = 0; 96 double kurtosis = 0; 97 for (int i = 0; i < u.size(); ++i) 98 { 99 double d = (u[i] - mean); 100 double d2 = sqr(d); 101 var += d2; 102 skew += d * 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 = std::exp(p.m() + sqr(p.s())/2); 111 double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s())); 112 double x_skew = (std::exp(sqr(p.s())) + 2) * 113 std::sqrt((std::exp(sqr(p.s())) - 1)); 114 double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) + 115 3*std::exp(2*sqr(p.s())) - 6; 116 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 117 assert(std::abs((var - x_var) / x_var) < 0.01); 118 assert(std::abs((skew - x_skew) / x_skew) < 0.01); 119 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03); 120 } 121 { 122 typedef std::lognormal_distribution<> D; 123 typedef D::param_type P; 124 typedef std::mt19937 G; 125 G g; 126 D d; 127 P p(-1./8, 0.5); 128 const int N = 1000000; 129 std::vector<D::result_type> u; 130 for (int i = 0; i < N; ++i) 131 { 132 D::result_type v = d(g, p); 133 assert(v > 0); 134 u.push_back(v); 135 } 136 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size(); 137 double var = 0; 138 double skew = 0; 139 double kurtosis = 0; 140 for (int i = 0; i < u.size(); ++i) 141 { 142 double d = (u[i] - mean); 143 double d2 = sqr(d); 144 var += d2; 145 skew += d * 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 = std::exp(p.m() + sqr(p.s())/2); 154 double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s())); 155 double x_skew = (std::exp(sqr(p.s())) + 2) * 156 std::sqrt((std::exp(sqr(p.s())) - 1)); 157 double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) + 158 3*std::exp(2*sqr(p.s())) - 6; 159 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 160 assert(std::abs((var - x_var) / x_var) < 0.01); 161 assert(std::abs((skew - x_skew) / x_skew) < 0.02); 162 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.05); 163 } 164 { 165 typedef std::lognormal_distribution<> D; 166 typedef D::param_type P; 167 typedef std::mt19937 G; 168 G g; 169 D d(3, 4); 170 P p; 171 const int N = 1000000; 172 std::vector<D::result_type> u; 173 for (int i = 0; i < N; ++i) 174 { 175 D::result_type v = d(g, p); 176 assert(v > 0); 177 u.push_back(v); 178 } 179 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size(); 180 double var = 0; 181 double skew = 0; 182 double kurtosis = 0; 183 for (int i = 0; i < u.size(); ++i) 184 { 185 double d = (u[i] - mean); 186 double d2 = sqr(d); 187 var += d2; 188 skew += d * d2; 189 kurtosis += d2 * d2; 190 } 191 var /= u.size(); 192 double dev = std::sqrt(var); 193 skew /= u.size() * dev * var; 194 kurtosis /= u.size() * var * var; 195 kurtosis -= 3; 196 double x_mean = std::exp(p.m() + sqr(p.s())/2); 197 double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s())); 198 double x_skew = (std::exp(sqr(p.s())) + 2) * 199 std::sqrt((std::exp(sqr(p.s())) - 1)); 200 double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) + 201 3*std::exp(2*sqr(p.s())) - 6; 202 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 203 assert(std::abs((var - x_var) / x_var) < 0.02); 204 assert(std::abs((skew - x_skew) / x_skew) < 0.08); 205 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.4); 206 } 207 { 208 typedef std::lognormal_distribution<> D; 209 typedef D::param_type P; 210 typedef std::mt19937 G; 211 G g; 212 D d; 213 P p(-0.78125, 1.25); 214 const int N = 1000000; 215 std::vector<D::result_type> u; 216 for (int i = 0; i < N; ++i) 217 { 218 D::result_type v = d(g, p); 219 assert(v > 0); 220 u.push_back(v); 221 } 222 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size(); 223 double var = 0; 224 double skew = 0; 225 double kurtosis = 0; 226 for (int i = 0; i < u.size(); ++i) 227 { 228 double d = (u[i] - mean); 229 double d2 = sqr(d); 230 var += d2; 231 skew += d * d2; 232 kurtosis += d2 * d2; 233 } 234 var /= u.size(); 235 double dev = std::sqrt(var); 236 skew /= u.size() * dev * var; 237 kurtosis /= u.size() * var * var; 238 kurtosis -= 3; 239 double x_mean = std::exp(p.m() + sqr(p.s())/2); 240 double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s())); 241 double x_skew = (std::exp(sqr(p.s())) + 2) * 242 std::sqrt((std::exp(sqr(p.s())) - 1)); 243 double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) + 244 3*std::exp(2*sqr(p.s())) - 6; 245 assert(std::abs((mean - x_mean) / x_mean) < 0.01); 246 assert(std::abs((var - x_var) / x_var) < 0.04); 247 assert(std::abs((skew - x_skew) / x_skew) < 0.2); 248 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.7); 249 } 250} 251