Home | History | Annotate | Download | only in rand.dist.pois.extreme
      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);
     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 std::mt19937 G;
     37     G g;
     38     D d(0.5, 2);
     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         u.push_back(v);
     45     }
     46     double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
     47     double var = 0;
     48     double skew = 0;
     49     double kurtosis = 0;
     50     for (unsigned i = 0; i < u.size(); ++i)
     51     {
     52         double dbl = (u[i] - mean);
     53         double d2 = sqr(dbl);
     54         var += d2;
     55         skew += dbl * d2;
     56         kurtosis += d2 * d2;
     57     }
     58     var /= u.size();
     59     double dev = std::sqrt(var);
     60     skew /= u.size() * dev * var;
     61     kurtosis /= u.size() * var * var;
     62     kurtosis -= 3;
     63     double x_mean = d.a() + d.b() * 0.577215665;
     64     double x_var = sqr(d.b()) * 1.644934067;
     65     double x_skew = 1.139547;
     66     double x_kurtosis = 12./5;
     67     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
     68     assert(std::abs((var - x_var) / x_var) < 0.01);
     69     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
     70     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
     71 }
     72 
     73 void
     74 test2()
     75 {
     76     typedef std::extreme_value_distribution<> D;
     77     typedef std::mt19937 G;
     78     G g;
     79     D d(1, 2);
     80     const int N = 1000000;
     81     std::vector<D::result_type> u;
     82     for (int i = 0; i < N; ++i)
     83     {
     84         D::result_type v = d(g);
     85         u.push_back(v);
     86     }
     87     double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
     88     double var = 0;
     89     double skew = 0;
     90     double kurtosis = 0;
     91     for (unsigned i = 0; i < u.size(); ++i)
     92     {
     93         double dbl = (u[i] - mean);
     94         double d2 = sqr(dbl);
     95         var += d2;
     96         skew += dbl * d2;
     97         kurtosis += d2 * d2;
     98     }
     99     var /= u.size();
    100     double dev = std::sqrt(var);
    101     skew /= u.size() * dev * var;
    102     kurtosis /= u.size() * var * var;
    103     kurtosis -= 3;
    104     double x_mean = d.a() + d.b() * 0.577215665;
    105     double x_var = sqr(d.b()) * 1.644934067;
    106     double x_skew = 1.139547;
    107     double x_kurtosis = 12./5;
    108     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
    109     assert(std::abs((var - x_var) / x_var) < 0.01);
    110     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
    111     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
    112 }
    113 
    114 void
    115 test3()
    116 {
    117     typedef std::extreme_value_distribution<> D;
    118     typedef std::mt19937 G;
    119     G g;
    120     D d(1.5, 3);
    121     const int N = 1000000;
    122     std::vector<D::result_type> u;
    123     for (int i = 0; i < N; ++i)
    124     {
    125         D::result_type v = d(g);
    126         u.push_back(v);
    127     }
    128     double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
    129     double var = 0;
    130     double skew = 0;
    131     double kurtosis = 0;
    132     for (unsigned i = 0; i < u.size(); ++i)
    133     {
    134         double dbl = (u[i] - mean);
    135         double d2 = sqr(dbl);
    136         var += d2;
    137         skew += dbl * d2;
    138         kurtosis += d2 * d2;
    139     }
    140     var /= u.size();
    141     double dev = std::sqrt(var);
    142     skew /= u.size() * dev * var;
    143     kurtosis /= u.size() * var * var;
    144     kurtosis -= 3;
    145     double x_mean = d.a() + d.b() * 0.577215665;
    146     double x_var = sqr(d.b()) * 1.644934067;
    147     double x_skew = 1.139547;
    148     double x_kurtosis = 12./5;
    149     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
    150     assert(std::abs((var - x_var) / x_var) < 0.01);
    151     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
    152     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
    153 }
    154 
    155 void
    156 test4()
    157 {
    158     typedef std::extreme_value_distribution<> D;
    159     typedef std::mt19937 G;
    160     G g;
    161     D d(3, 4);
    162     const int N = 1000000;
    163     std::vector<D::result_type> u;
    164     for (int i = 0; i < N; ++i)
    165     {
    166         D::result_type v = d(g);
    167         u.push_back(v);
    168     }
    169     double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
    170     double var = 0;
    171     double skew = 0;
    172     double kurtosis = 0;
    173     for (unsigned i = 0; i < u.size(); ++i)
    174     {
    175         double dbl = (u[i] - mean);
    176         double d2 = sqr(dbl);
    177         var += d2;
    178         skew += dbl * d2;
    179         kurtosis += d2 * d2;
    180     }
    181     var /= u.size();
    182     double dev = std::sqrt(var);
    183     skew /= u.size() * dev * var;
    184     kurtosis /= u.size() * var * var;
    185     kurtosis -= 3;
    186     double x_mean = d.a() + d.b() * 0.577215665;
    187     double x_var = sqr(d.b()) * 1.644934067;
    188     double x_skew = 1.139547;
    189     double x_kurtosis = 12./5;
    190     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
    191     assert(std::abs((var - x_var) / x_var) < 0.01);
    192     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
    193     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
    194 }
    195 
    196 int main()
    197 {
    198     test1();
    199     test2();
    200     test3();
    201     test4();
    202 }
    203