test/std/numerics/rand/rand.dis/rand.dist.pois/rand.dist.pois.extreme/eval_param.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 RealType = double>
  15 // class extreme_value_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
  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::extreme_value_distribution<> D;
  36     typedef D::param_type P;
  37     typedef std::mt19937 G;
  38     G g;
  39     D d(-0.5, 1);
  40     P p(0.5, 2);
  41     const int N = 1000000;
  42     std::vector<D::result_type> u;
  43     for (int i = 0; i < N; ++i)
  44     {
  45         D::result_type v = d(g, p);
  46         u.push_back(v);
  47     }
  48     double mean = std::accumulate(u.begin(), u.end(), 0.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 = p.a() + p.b() * 0.577215665;
  66     double x_var = sqr(p.b()) * 1.644934067;
  67     double x_skew = 1.139547;
  68     double x_kurtosis = 12./5;
  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.01);
  73 }
  74
  75 void
  76 test2()
  77 {
  78     typedef std::extreme_value_distribution<> D;
  79     typedef D::param_type P;
  80     typedef std::mt19937 G;
  81     G g;
  82     D d(-0.5, 1);
  83     P p(1, 2);
  84     const int N = 1000000;
  85     std::vector<D::result_type> u;
  86     for (int i = 0; i < N; ++i)
  87     {
  88         D::result_type v = d(g, p);
  89         u.push_back(v);
  90     }
  91     double mean = std::accumulate(u.begin(), u.end(), 0.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 = p.a() + p.b() * 0.577215665;
 109     double x_var = sqr(p.b()) * 1.644934067;
 110     double x_skew = 1.139547;
 111     double x_kurtosis = 12./5;
 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::extreme_value_distribution<> D;
 122     typedef D::param_type P;
 123     typedef std::mt19937 G;
 124     G g;
 125     D d(-0.5, 1);
 126     P p(1.5, 3);
 127     const int N = 1000000;
 128     std::vector<D::result_type> u;
 129     for (int i = 0; i < N; ++i)
 130     {
 131         D::result_type v = d(g, p);
 132         u.push_back(v);
 133     }
 134     double mean = std::accumulate(u.begin(), u.end(), 0.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 = p.a() + p.b() * 0.577215665;
 152     double x_var = sqr(p.b()) * 1.644934067;
 153     double x_skew = 1.139547;
 154     double x_kurtosis = 12./5;
 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.01);
 159 }
 160
 161 void
 162 test4()
 163 {
 164     typedef std::extreme_value_distribution<> D;
 165     typedef D::param_type P;
 166     typedef std::mt19937 G;
 167     G g;
 168     D d(-0.5, 1);
 169     P p(3, 4);
 170     const int N = 1000000;
 171     std::vector<D::result_type> u;
 172     for (int i = 0; i < N; ++i)
 173     {
 174         D::result_type v = d(g, p);
 175         u.push_back(v);
 176     }
 177     double mean = std::accumulate(u.begin(), u.end(), 0.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 = p.a() + p.b() * 0.577215665;
 195     double x_var = sqr(p.b()) * 1.644934067;
 196     double x_skew = 1.139547;
 197     double x_kurtosis = 12./5;
 198     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 199     assert(std::abs((var - x_var) / x_var) < 0.01);
 200     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
 201     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
 202 }
 203
 204 int main()
 205 {
 206     test1();
 207     test2();
 208     test3();
 209     test4();
 210 }