test/std/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.negbin/eval.pass.cpp

   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 IntType = int>
  15 // class negative_binomial_distribution
  16
  17 // template<class _URNG> result_type operator()(_URNG& g);
  18
  19 #include <random>
  20 #include <numeric>
  21 #include <vector>
  22 #include <cassert>
  23
  24 template <class T>
  25 inline
  26 T
  27 sqr(T x)
  28 {
  29     return x * x;
  30 }
  31
  32 void
  33 test1()
  34 {
  35     typedef std::negative_binomial_distribution<> D;
  36     typedef std::minstd_rand G;
  37     G g;
  38     D d(5, .25);
  39     const int N = 1000000;
  40     std::vector<D::result_type> u;
  41     for (int i = 0; i < N; ++i)
  42     {
  43         D::result_type v = d(g);
  44         assert(d.min() <= v && v <= d.max());
  45         u.push_back(v);
  46     }
  47     double mean = std::accumulate(u.begin(), u.end(),
  48                                           double(0)) / u.size();
  49     double var = 0;
  50     double skew = 0;
  51     double kurtosis = 0;
  52     for (unsigned i = 0; i < u.size(); ++i)
  53     {
  54         double dbl = (u[i] - mean);
  55         double d2 = sqr(dbl);
  56         var += d2;
  57         skew += dbl * d2;
  58         kurtosis += d2 * d2;
  59     }
  60     var /= u.size();
  61     double dev = std::sqrt(var);
  62     skew /= u.size() * dev * var;
  63     kurtosis /= u.size() * var * var;
  64     kurtosis -= 3;
  65     double x_mean = d.k() * (1 - d.p()) / d.p();
  66     double x_var = x_mean / d.p();
  67     double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
  68     double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
  69     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
  70     assert(std::abs((var - x_var) / x_var) < 0.01);
  71     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
  72     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
  73 }
  74
  75 void
  76 test2()
  77 {
  78     typedef std::negative_binomial_distribution<> D;
  79     typedef std::mt19937 G;
  80     G g;
  81     D d(30, .03125);
  82     const int N = 1000000;
  83     std::vector<D::result_type> u;
  84     for (int i = 0; i < N; ++i)
  85     {
  86         D::result_type v = d(g);
  87         assert(d.min() <= v && v <= d.max());
  88         u.push_back(v);
  89     }
  90     double mean = std::accumulate(u.begin(), u.end(),
  91                                           double(0)) / u.size();
  92     double var = 0;
  93     double skew = 0;
  94     double kurtosis = 0;
  95     for (unsigned i = 0; i < u.size(); ++i)
  96     {
  97         double dbl = (u[i] - mean);
  98         double d2 = sqr(dbl);
  99         var += d2;
 100         skew += dbl * d2;
 101         kurtosis += d2 * d2;
 102     }
 103     var /= u.size();
 104     double dev = std::sqrt(var);
 105     skew /= u.size() * dev * var;
 106     kurtosis /= u.size() * var * var;
 107     kurtosis -= 3;
 108     double x_mean = d.k() * (1 - d.p()) / d.p();
 109     double x_var = x_mean / d.p();
 110     double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
 111     double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
 112     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 113     assert(std::abs((var - x_var) / x_var) < 0.01);
 114     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
 115     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
 116 }
 117
 118 void
 119 test3()
 120 {
 121     typedef std::negative_binomial_distribution<> D;
 122     typedef std::mt19937 G;
 123     G g;
 124     D d(40, .25);
 125     const int N = 1000000;
 126     std::vector<D::result_type> u;
 127     for (int i = 0; i < N; ++i)
 128     {
 129         D::result_type v = d(g);
 130         assert(d.min() <= v && v <= d.max());
 131         u.push_back(v);
 132     }
 133     double mean = std::accumulate(u.begin(), u.end(),
 134                                           double(0)) / u.size();
 135     double var = 0;
 136     double skew = 0;
 137     double kurtosis = 0;
 138     for (unsigned i = 0; i < u.size(); ++i)
 139     {
 140         double dbl = (u[i] - mean);
 141         double d2 = sqr(dbl);
 142         var += d2;
 143         skew += dbl * d2;
 144         kurtosis += d2 * d2;
 145     }
 146     var /= u.size();
 147     double dev = std::sqrt(var);
 148     skew /= u.size() * dev * var;
 149     kurtosis /= u.size() * var * var;
 150     kurtosis -= 3;
 151     double x_mean = d.k() * (1 - d.p()) / d.p();
 152     double x_var = x_mean / d.p();
 153     double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
 154     double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
 155     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 156     assert(std::abs((var - x_var) / x_var) < 0.01);
 157     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
 158     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
 159 }
 160
 161 void
 162 test4()
 163 {
 164     typedef std::negative_binomial_distribution<> D;
 165     typedef std::mt19937 G;
 166     G g;
 167     D d(40, 1);
 168     const int N = 1000;
 169     std::vector<D::result_type> u;
 170     for (int i = 0; i < N; ++i)
 171     {
 172         D::result_type v = d(g);
 173         assert(d.min() <= v && v <= d.max());
 174         u.push_back(v);
 175     }
 176     double mean = std::accumulate(u.begin(), u.end(),
 177                                           double(0)) / u.size();
 178     double var = 0;
 179     double skew = 0;
 180     double kurtosis = 0;
 181     for (unsigned i = 0; i < u.size(); ++i)
 182     {
 183         double dbl = (u[i] - mean);
 184         double d2 = sqr(dbl);
 185         var += d2;
 186         skew += dbl * d2;
 187         kurtosis += d2 * d2;
 188     }
 189     var /= u.size();
 190     double dev = std::sqrt(var);
 191     skew /= u.size() * dev * var;
 192     kurtosis /= u.size() * var * var;
 193     kurtosis -= 3;
 194     double x_mean = d.k() * (1 - d.p()) / d.p();
 195     double x_var = x_mean / d.p();
 196 //    double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
 197 //    double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
 198     assert(mean == x_mean);
 199     assert(var == x_var);
 200 }
 201
 202 void
 203 test5()
 204 {
 205     typedef std::negative_binomial_distribution<> D;
 206     typedef std::mt19937 G;
 207     G g;
 208     D d(400, 0.5);
 209     const int N = 1000000;
 210     std::vector<D::result_type> u;
 211     for (int i = 0; i < N; ++i)
 212     {
 213         D::result_type v = d(g);
 214         assert(d.min() <= v && v <= d.max());
 215         u.push_back(v);
 216     }
 217     double mean = std::accumulate(u.begin(), u.end(),
 218                                           double(0)) / u.size();
 219     double var = 0;
 220     double skew = 0;
 221     double kurtosis = 0;
 222     for (unsigned i = 0; i < u.size(); ++i)
 223     {
 224         double dbl = (u[i] - mean);
 225         double d2 = sqr(dbl);
 226         var += d2;
 227         skew += dbl * d2;
 228         kurtosis += d2 * d2;
 229     }
 230     var /= u.size();
 231     double dev = std::sqrt(var);
 232     skew /= u.size() * dev * var;
 233     kurtosis /= u.size() * var * var;
 234     kurtosis -= 3;
 235     double x_mean = d.k() * (1 - d.p()) / d.p();
 236     double x_var = x_mean / d.p();
 237     double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
 238     double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
 239     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 240     assert(std::abs((var - x_var) / x_var) < 0.01);
 241     assert(std::abs((skew - x_skew) / x_skew) < 0.04);
 242     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.05);
 243 }
 244
 245 void
 246 test6()
 247 {
 248     typedef std::negative_binomial_distribution<> D;
 249     typedef std::mt19937 G;
 250     G g;
 251     D d(1, 0.05);
 252     const int N = 1000000;
 253     std::vector<D::result_type> u;
 254     for (int i = 0; i < N; ++i)
 255     {
 256         D::result_type v = d(g);
 257         assert(d.min() <= v && v <= d.max());
 258         u.push_back(v);
 259     }
 260     double mean = std::accumulate(u.begin(), u.end(),
 261                                           double(0)) / u.size();
 262     double var = 0;
 263     double skew = 0;
 264     double kurtosis = 0;
 265     for (unsigned i = 0; i < u.size(); ++i)
 266     {
 267         double dbl = (u[i] - mean);
 268         double d2 = sqr(dbl);
 269         var += d2;
 270         skew += dbl * d2;
 271         kurtosis += d2 * d2;
 272     }
 273     var /= u.size();
 274     double dev = std::sqrt(var);
 275     skew /= u.size() * dev * var;
 276     kurtosis /= u.size() * var * var;
 277     kurtosis -= 3;
 278     double x_mean = d.k() * (1 - d.p()) / d.p();
 279     double x_var = x_mean / d.p();
 280     double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
 281     double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
 282     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 283     assert(std::abs((var - x_var) / x_var) < 0.01);
 284     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
 285     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
 286 }
 287
 288 int main()
 289 {
 290     test1();
 291     test2();
 292     test3();
 293     test4();
 294     test5();
 295     test6();
 296 }