libcxx/test/std/numerics/rand/rand.dist/rand.dist.pois/rand.dist.pois.extreme/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 // <random>
  12
  13 // template<class RealType = double>
  14 // class extreme_value_distribution
  15
  16 // template<class _URNG> result_type operator()(_URNG& g);
  17
  18 #include <random>
  19 #include <cassert>
  20 #include <cmath>
  21 #include <numeric>
  22 #include <vector>
  23
  24 #include "test_macros.h"
  25
  26 template <class T>
  27 inline
  28 T
  29 sqr(T x)
  30 {
  31     return x * x;
  32 }
  33
  34 void
  35 test1()
  36 {
  37     typedef std::extreme_value_distribution<> D;
  38     typedef std::mt19937 G;
  39     G g;
  40     D d(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);
  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 = d.a() + d.b() * 0.577215665;
  66     double x_var = sqr(d.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.03);
  73 }
  74
  75 void
  76 test2()
  77 {
  78     typedef std::extreme_value_distribution<> D;
  79     typedef std::mt19937 G;
  80     G g;
  81     D d(1, 2);
  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         u.push_back(v);
  88     }
  89     double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
  90     double var = 0;
  91     double skew = 0;
  92     double kurtosis = 0;
  93     for (unsigned i = 0; i < u.size(); ++i)
  94     {
  95         double dbl = (u[i] - mean);
  96         double d2 = sqr(dbl);
  97         var += d2;
  98         skew += dbl * d2;
  99         kurtosis += d2 * d2;
 100     }
 101     var /= u.size();
 102     double dev = std::sqrt(var);
 103     skew /= u.size() * dev * var;
 104     kurtosis /= u.size() * var * var;
 105     kurtosis -= 3;
 106     double x_mean = d.a() + d.b() * 0.577215665;
 107     double x_var = sqr(d.b()) * 1.644934067;
 108     double x_skew = 1.139547;
 109     double x_kurtosis = 12./5;
 110     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 111     assert(std::abs((var - x_var) / x_var) < 0.01);
 112     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
 113     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
 114 }
 115
 116 void
 117 test3()
 118 {
 119     typedef std::extreme_value_distribution<> D;
 120     typedef std::mt19937 G;
 121     G g;
 122     D d(1.5, 3);
 123     const int N = 1000000;
 124     std::vector<D::result_type> u;
 125     for (int i = 0; i < N; ++i)
 126     {
 127         D::result_type v = d(g);
 128         u.push_back(v);
 129     }
 130     double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
 131     double var = 0;
 132     double skew = 0;
 133     double kurtosis = 0;
 134     for (unsigned i = 0; i < u.size(); ++i)
 135     {
 136         double dbl = (u[i] - mean);
 137         double d2 = sqr(dbl);
 138         var += d2;
 139         skew += dbl * d2;
 140         kurtosis += d2 * d2;
 141     }
 142     var /= u.size();
 143     double dev = std::sqrt(var);
 144     skew /= u.size() * dev * var;
 145     kurtosis /= u.size() * var * var;
 146     kurtosis -= 3;
 147     double x_mean = d.a() + d.b() * 0.577215665;
 148     double x_var = sqr(d.b()) * 1.644934067;
 149     double x_skew = 1.139547;
 150     double x_kurtosis = 12./5;
 151     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 152     assert(std::abs((var - x_var) / x_var) < 0.01);
 153     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
 154     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
 155 }
 156
 157 void
 158 test4()
 159 {
 160     typedef std::extreme_value_distribution<> D;
 161     typedef std::mt19937 G;
 162     G g;
 163     D d(3, 4);
 164     const int N = 1000000;
 165     std::vector<D::result_type> u;
 166     for (int i = 0; i < N; ++i)
 167     {
 168         D::result_type v = d(g);
 169         u.push_back(v);
 170     }
 171     double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
 172     double var = 0;
 173     double skew = 0;
 174     double kurtosis = 0;
 175     for (unsigned i = 0; i < u.size(); ++i)
 176     {
 177         double dbl = (u[i] - mean);
 178         double d2 = sqr(dbl);
 179         var += d2;
 180         skew += dbl * d2;
 181         kurtosis += d2 * d2;
 182     }
 183     var /= u.size();
 184     double dev = std::sqrt(var);
 185     skew /= u.size() * dev * var;
 186     kurtosis /= u.size() * var * var;
 187     kurtosis -= 3;
 188     double x_mean = d.a() + d.b() * 0.577215665;
 189     double x_var = sqr(d.b()) * 1.644934067;
 190     double x_skew = 1.139547;
 191     double x_kurtosis = 12./5;
 192     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
 193     assert(std::abs((var - x_var) / x_var) < 0.01);
 194     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
 195     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
 196 }
 197
 198 int main(int, char**)
 199 {
 200     test1();
 201     test2();
 202     test3();
 203     test4();
 204
 205   return 0;
 206 }