Home | History | Annotate | Download | only in rand.dist.pois.exp
      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 exponential_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::exponential_distribution<> D;
     36         typedef D::param_type P;
     37         typedef std::mt19937 G;
     38         G g;
     39         D d(.75);
     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             assert(d.min() < v);
     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 (int i = 0; i < u.size(); ++i)
     53         {
     54             double d = (u[i] - mean);
     55             double d2 = sqr(d);
     56             var += d2;
     57             skew += d * 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 = 1/d.lambda();
     66         double x_var = 1/sqr(d.lambda());
     67         double x_skew = 2;
     68         double x_kurtosis = 6;
     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         typedef std::exponential_distribution<> D;
     76         typedef D::param_type P;
     77         typedef std::mt19937 G;
     78         G g;
     79         D d(1);
     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             assert(d.min() < v);
     86             u.push_back(v);
     87         }
     88         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
     89         double var = 0;
     90         double skew = 0;
     91         double kurtosis = 0;
     92         for (int i = 0; i < u.size(); ++i)
     93         {
     94             double d = (u[i] - mean);
     95             double d2 = sqr(d);
     96             var += d2;
     97             skew += d * d2;
     98             kurtosis += d2 * d2;
     99         }
    100         var /= u.size();
    101         double dev = std::sqrt(var);
    102         skew /= u.size() * dev * var;
    103         kurtosis /= u.size() * var * var;
    104         kurtosis -= 3;
    105         double x_mean = 1/d.lambda();
    106         double x_var = 1/sqr(d.lambda());
    107         double x_skew = 2;
    108         double x_kurtosis = 6;
    109         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
    110         assert(std::abs((var - x_var) / x_var) < 0.01);
    111         assert(std::abs((skew - x_skew) / x_skew) < 0.01);
    112         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
    113     }
    114     {
    115         typedef std::exponential_distribution<> D;
    116         typedef D::param_type P;
    117         typedef std::mt19937 G;
    118         G g;
    119         D d(10);
    120         const int N = 1000000;
    121         std::vector<D::result_type> u;
    122         for (int i = 0; i < N; ++i)
    123         {
    124             D::result_type v = d(g);
    125             assert(d.min() < v);
    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 (int i = 0; i < u.size(); ++i)
    133         {
    134             double d = (u[i] - mean);
    135             double d2 = sqr(d);
    136             var += d2;
    137             skew += d * 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 = 1/d.lambda();
    146         double x_var = 1/sqr(d.lambda());
    147         double x_skew = 2;
    148         double x_kurtosis = 6;
    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