libcxx/test/std/numerics/rand/rand.dist/rand.dist.bern/rand.dist.bern.negbin/eval.pass.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-exception
   6 //
   7 //===----------------------------------------------------------------------===//
   8 //
   9 // REQUIRES: long_tests
  10
  11 // This test is super slow, in particular with msan or tsan. In order to avoid timeouts and to
  12 // spend less time waiting for this particular test to complete we compile with optimizations.
  13 // ADDITIONAL_COMPILE_FLAGS(msan): -O1
  14 // ADDITIONAL_COMPILE_FLAGS(tsan): -O1
  15
  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_distribution
  23
  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
  34 template <class T>
  35 T sqr(T x) {
  36     return x * x;
  37 }
  38
  39 template <class T>
  40 void 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
  81 template <class T>
  82 void 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
 123 template <class T>
 124 void 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
 165 template <class T>
 166 void 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
 209 template <class T>
 210 void 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
 251 template <class T>
 252 void 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
 293 template <class T>
 294 void tests() {
 295     test1<T>();
 296     test2<T>();
 297     test3<T>();
 298     test4<T>();
 299     test5<T>();
 300     test6<T>();
 301 }
 302
 303 int 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) // extension
 315     // 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 #endif
 322 #endif
 323
 324     return 0;
 325 }