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