326 lines · cpp
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-exception6//7//===----------------------------------------------------------------------===//8//9// REQUIRES: long_tests10 11// This test is super slow, in particular with msan or tsan. In order to avoid timeouts and to12// spend less time waiting for this particular test to complete we compile with optimizations.13// ADDITIONAL_COMPILE_FLAGS(msan): -O114// ADDITIONAL_COMPILE_FLAGS(tsan): -O115 16// FIXME: This and other tests fail under GCC with optimizations enabled.17// More investigation is needed, but it appears that GCC is performing more constant folding.18 19// <random>20 21// template<class IntType = int>22// class negative_binomial_distribution23 24// template<class _URNG> result_type operator()(_URNG& g);25 26#include <random>27#include <cassert>28#include <cmath>29#include <numeric>30#include <vector>31 32#include "test_macros.h"33 34template <class T>35T sqr(T x) {36 return x * x;37}38 39template <class T>40void test1() {41 typedef std::negative_binomial_distribution<T> D;42 typedef std::minstd_rand G;43 G g;44 D d(5, .25);45 const int N = 1000000;46 std::vector<typename D::result_type> u;47 for (int i = 0; i < N; ++i)48 {49 typename D::result_type v = d(g);50 assert(d.min() <= v && v <= d.max());51 u.push_back(v);52 }53 double mean = std::accumulate(u.begin(), u.end(),54 double(0)) / u.size();55 double var = 0;56 double skew = 0;57 double kurtosis = 0;58 for (unsigned i = 0; i < u.size(); ++i)59 {60 double dbl = (u[i] - mean);61 double d2 = sqr(dbl);62 var += d2;63 skew += dbl * d2;64 kurtosis += d2 * d2;65 }66 var /= u.size();67 double dev = std::sqrt(var);68 skew /= u.size() * dev * var;69 kurtosis /= u.size() * var * var;70 kurtosis -= 3;71 double x_mean = d.k() * (1 - d.p()) / d.p();72 double x_var = x_mean / d.p();73 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));74 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));75 assert(std::abs((mean - x_mean) / x_mean) < 0.01);76 assert(std::abs((var - x_var) / x_var) < 0.01);77 assert(std::abs((skew - x_skew) / x_skew) < 0.01);78 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);79}80 81template <class T>82void test2() {83 typedef std::negative_binomial_distribution<T> D;84 typedef std::mt19937 G;85 G g;86 D d(30, .03125);87 const int N = 1000000;88 std::vector<typename D::result_type> u;89 for (int i = 0; i < N; ++i)90 {91 typename D::result_type v = d(g);92 assert(d.min() <= v && v <= d.max());93 u.push_back(v);94 }95 double mean = std::accumulate(u.begin(), u.end(),96 double(0)) / u.size();97 double var = 0;98 double skew = 0;99 double kurtosis = 0;100 for (unsigned i = 0; i < u.size(); ++i)101 {102 double dbl = (u[i] - mean);103 double d2 = sqr(dbl);104 var += d2;105 skew += dbl * d2;106 kurtosis += d2 * d2;107 }108 var /= u.size();109 double dev = std::sqrt(var);110 skew /= u.size() * dev * var;111 kurtosis /= u.size() * var * var;112 kurtosis -= 3;113 double x_mean = d.k() * (1 - d.p()) / d.p();114 double x_var = x_mean / d.p();115 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));116 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));117 assert(std::abs((mean - x_mean) / x_mean) < 0.01);118 assert(std::abs((var - x_var) / x_var) < 0.01);119 assert(std::abs((skew - x_skew) / x_skew) < 0.02);120 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.1);121}122 123template <class T>124void test3() {125 typedef std::negative_binomial_distribution<T> D;126 typedef std::mt19937 G;127 G g;128 D d(40, .25);129 const int N = 1000000;130 std::vector<typename D::result_type> u;131 for (int i = 0; i < N; ++i)132 {133 typename D::result_type v = d(g);134 assert(d.min() <= v && v <= d.max());135 u.push_back(v);136 }137 double mean = std::accumulate(u.begin(), u.end(),138 double(0)) / u.size();139 double var = 0;140 double skew = 0;141 double kurtosis = 0;142 for (unsigned i = 0; i < u.size(); ++i)143 {144 double dbl = (u[i] - mean);145 double d2 = sqr(dbl);146 var += d2;147 skew += dbl * d2;148 kurtosis += d2 * d2;149 }150 var /= u.size();151 double dev = std::sqrt(var);152 skew /= u.size() * dev * var;153 kurtosis /= u.size() * var * var;154 kurtosis -= 3;155 double x_mean = d.k() * (1 - d.p()) / d.p();156 double x_var = x_mean / d.p();157 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));158 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));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.08);163}164 165template <class T>166void test4() {167 typedef std::negative_binomial_distribution<T> D;168 typedef std::mt19937 G;169 G g;170 D d(40, 1);171 const int N = 1000;172 std::vector<typename D::result_type> u;173 for (int i = 0; i < N; ++i)174 {175 typename D::result_type v = d(g);176 assert(d.min() <= v && v <= d.max());177 u.push_back(v);178 }179 double mean = std::accumulate(u.begin(), u.end(),180 double(0)) / u.size();181 double var = 0;182 double skew = 0;183 double kurtosis = 0;184 for (unsigned i = 0; i < u.size(); ++i)185 {186 double dbl = (u[i] - mean);187 double d2 = sqr(dbl);188 var += d2;189 skew += dbl * d2;190 kurtosis += d2 * d2;191 }192 var /= u.size();193 double dev = std::sqrt(var);194 skew /= u.size() * dev * var;195 kurtosis /= u.size() * var * var;196 kurtosis -= 3;197 double x_mean = d.k() * (1 - d.p()) / d.p();198 double x_var = x_mean / d.p();199 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));200 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));201 assert(mean == x_mean);202 assert(var == x_var);203 // assert(skew == x_skew);204 (void)skew; (void)x_skew;205 // assert(kurtosis == x_kurtosis);206 (void)kurtosis; (void)x_kurtosis;207}208 209template <class T>210void test5() {211 typedef std::negative_binomial_distribution<T> D;212 typedef std::mt19937 G;213 G g;214 D d(127, 0.5);215 const int N = 1000000;216 std::vector<typename D::result_type> u;217 for (int i = 0; i < N; ++i)218 {219 typename D::result_type v = d(g);220 assert(d.min() <= v && v <= d.max());221 u.push_back(v);222 }223 double mean = std::accumulate(u.begin(), u.end(),224 double(0)) / u.size();225 double var = 0;226 double skew = 0;227 double kurtosis = 0;228 for (unsigned i = 0; i < u.size(); ++i)229 {230 double dbl = (u[i] - mean);231 double d2 = sqr(dbl);232 var += d2;233 skew += dbl * d2;234 kurtosis += d2 * d2;235 }236 var /= u.size();237 double dev = std::sqrt(var);238 skew /= u.size() * dev * var;239 kurtosis /= u.size() * var * var;240 kurtosis -= 3;241 double x_mean = d.k() * (1 - d.p()) / d.p();242 double x_var = x_mean / d.p();243 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));244 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));245 assert(std::abs((mean - x_mean) / x_mean) < 0.01);246 assert(std::abs((var - x_var) / x_var) < 0.01);247 assert(std::abs((skew - x_skew) / x_skew) < 0.02);248 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.3);249}250 251template <class T>252void test6() {253 typedef std::negative_binomial_distribution<T> D;254 typedef std::mt19937 G;255 G g;256 D d(1, 0.05);257 const int N = 1000000;258 std::vector<typename D::result_type> u;259 for (int i = 0; i < N; ++i)260 {261 typename D::result_type v = d(g);262 assert(d.min() <= v && v <= d.max());263 u.push_back(v);264 }265 double mean = std::accumulate(u.begin(), u.end(),266 double(0)) / u.size();267 double var = 0;268 double skew = 0;269 double kurtosis = 0;270 for (unsigned i = 0; i < u.size(); ++i)271 {272 double dbl = (u[i] - mean);273 double d2 = sqr(dbl);274 var += d2;275 skew += dbl * d2;276 kurtosis += d2 * d2;277 }278 var /= u.size();279 double dev = std::sqrt(var);280 skew /= u.size() * dev * var;281 kurtosis /= u.size() * var * var;282 kurtosis -= 3;283 double x_mean = d.k() * (1 - d.p()) / d.p();284 double x_var = x_mean / d.p();285 double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));286 double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));287 assert(std::abs((mean - x_mean) / x_mean) < 0.01);288 assert(std::abs((var - x_var) / x_var) < 0.01);289 assert(std::abs((skew - x_skew) / x_skew) < 0.01);290 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);291}292 293template <class T>294void tests() {295 test1<T>();296 test2<T>();297 test3<T>();298 test4<T>();299 test5<T>();300 test6<T>();301}302 303int main(int, char**) {304 tests<short>();305 tests<int>();306 tests<long>();307 tests<long long>();308 309 tests<unsigned short>();310 tests<unsigned int>();311 tests<unsigned long>();312 tests<unsigned long long>();313 314#if defined(_LIBCPP_VERSION) // extension315 // TODO: std::negative_binomial_distribution currently doesn't work reliably with small types.316 // tests<int8_t>();317 // tests<uint8_t>();318#if !defined(TEST_HAS_NO_INT128)319 tests<__int128_t>();320 tests<__uint128_t>();321#endif322#endif323 324 return 0;325}326