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 #include <cstddef>
     24 
     25 template <class T>
     26 inline
     27 T
     28 sqr(T x)
     29 {
     30     return x * x;
     31 }
     32 
     33 int main()
     34 {
     35     {
     36         typedef std::exponential_distribution<> D;
     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 (std::size_t 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 = 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 std::mt19937 G;
     77         G g;
     78         D d(1);
     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             assert(d.min() < v);
     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 (std::size_t 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 = 1/d.lambda();
    105         double x_var = 1/sqr(d.lambda());
    106         double x_skew = 2;
    107         double x_kurtosis = 6;
    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         typedef std::exponential_distribution<> D;
    115         typedef std::mt19937 G;
    116         G g;
    117         D d(10);
    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             assert(d.min() < v);
    124             u.push_back(v);
    125         }
    126         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
    127         double var = 0;
    128         double skew = 0;
    129         double kurtosis = 0;
    130         for (std::size_t i = 0; i < u.size(); ++i)
    131         {
    132             double dbl = (u[i] - mean);
    133             double d2 = sqr(dbl);
    134             var += d2;
    135             skew += dbl * d2;
    136             kurtosis += d2 * d2;
    137         }
    138         var /= u.size();
    139         double dev = std::sqrt(var);
    140         skew /= u.size() * dev * var;
    141         kurtosis /= u.size() * var * var;
    142         kurtosis -= 3;
    143         double x_mean = 1/d.lambda();
    144         double x_var = 1/sqr(d.lambda());
    145         double x_skew = 2;
    146         double x_kurtosis = 6;
    147         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
    148         assert(std::abs((var - x_var) / x_var) < 0.01);
    149         assert(std::abs((skew - x_skew) / x_skew) < 0.01);
    150         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
    151     }
    152 }
    153