Home | History | Annotate | Download | only in rand.dist.bern.geo
      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 IntType = int>
     15 // class geometric_distribution
     16 
     17 // template<class _URNG> result_type operator()(_URNG& g);
     18 
     19 #include <random>
     20 #include <numeric>
     21 #include <vector>
     22 #include <cassert>
     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::geometric_distribution<> D;
     36         typedef std::mt19937 G;
     37         G g;
     38         D d(.03125);
     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             assert(d.min() <= v && v <= d.max());
     45             u.push_back(v);
     46         }
     47         double mean = std::accumulate(u.begin(), u.end(),
     48                                               double(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.p()) / d.p();
     66         double x_var = x_mean / d.p();
     67         double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
     68         double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
     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::geometric_distribution<> D;
     76         typedef std::mt19937 G;
     77         G g;
     78         D d(0.05);
     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 && v <= d.max());
     85             u.push_back(v);
     86         }
     87         double mean = std::accumulate(u.begin(), u.end(),
     88                                               double(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.p()) / d.p();
    106         double x_var = x_mean / d.p();
    107         double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
    108         double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
    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.03);
    113     }
    114     {
    115         typedef std::geometric_distribution<> D;
    116         typedef std::minstd_rand G;
    117         G g;
    118         D d(.25);
    119         const int N = 1000000;
    120         std::vector<D::result_type> u;
    121         for (int i = 0; i < N; ++i)
    122         {
    123             D::result_type v = d(g);
    124             assert(d.min() <= v && v <= d.max());
    125             u.push_back(v);
    126         }
    127         double mean = std::accumulate(u.begin(), u.end(),
    128                                               double(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.p()) / d.p();
    146         double x_var = x_mean / d.p();
    147         double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
    148         double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
    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.02);
    153     }
    154     {
    155         typedef std::geometric_distribution<> D;
    156         typedef std::mt19937 G;
    157         G g;
    158         D d(0.5);
    159         const int N = 1000000;
    160         std::vector<D::result_type> u;
    161         for (int i = 0; i < N; ++i)
    162         {
    163             D::result_type v = d(g);
    164             assert(d.min() <= v && v <= d.max());
    165             u.push_back(v);
    166         }
    167         double mean = std::accumulate(u.begin(), u.end(),
    168                                               double(0)) / u.size();
    169         double var = 0;
    170         double skew = 0;
    171         double kurtosis = 0;
    172         for (int i = 0; i < u.size(); ++i)
    173         {
    174             double d = (u[i] - mean);
    175             double d2 = sqr(d);
    176             var += d2;
    177             skew += d * d2;
    178             kurtosis += d2 * d2;
    179         }
    180         var /= u.size();
    181         double dev = std::sqrt(var);
    182         skew /= u.size() * dev * var;
    183         kurtosis /= u.size() * var * var;
    184         kurtosis -= 3;
    185         double x_mean = (1 - d.p()) / d.p();
    186         double x_var = x_mean / d.p();
    187         double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
    188         double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
    189         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
    190         assert(std::abs((var - x_var) / x_var) < 0.01);
    191         assert(std::abs((skew - x_skew) / x_skew) < 0.01);
    192         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
    193     }
    194     {
    195         typedef std::geometric_distribution<> D;
    196         typedef std::mt19937 G;
    197         G g;
    198         D d(0.75);
    199         const int N = 1000000;
    200         std::vector<D::result_type> u;
    201         for (int i = 0; i < N; ++i)
    202         {
    203             D::result_type v = d(g);
    204             assert(d.min() <= v && v <= d.max());
    205             u.push_back(v);
    206         }
    207         double mean = std::accumulate(u.begin(), u.end(),
    208                                               double(0)) / u.size();
    209         double var = 0;
    210         double skew = 0;
    211         double kurtosis = 0;
    212         for (int i = 0; i < u.size(); ++i)
    213         {
    214             double d = (u[i] - mean);
    215             double d2 = sqr(d);
    216             var += d2;
    217             skew += d * d2;
    218             kurtosis += d2 * d2;
    219         }
    220         var /= u.size();
    221         double dev = std::sqrt(var);
    222         skew /= u.size() * dev * var;
    223         kurtosis /= u.size() * var * var;
    224         kurtosis -= 3;
    225         double x_mean = (1 - d.p()) / d.p();
    226         double x_var = x_mean / d.p();
    227         double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
    228         double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
    229         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
    230         assert(std::abs((var - x_var) / x_var) < 0.01);
    231         assert(std::abs((skew - x_skew) / x_skew) < 0.01);
    232         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
    233     }
    234     {
    235         typedef std::geometric_distribution<> D;
    236         typedef std::mt19937 G;
    237         G g;
    238         D d(0.96875);
    239         const int N = 1000000;
    240         std::vector<D::result_type> u;
    241         for (int i = 0; i < N; ++i)
    242         {
    243             D::result_type v = d(g);
    244             assert(d.min() <= v && v <= d.max());
    245             u.push_back(v);
    246         }
    247         double mean = std::accumulate(u.begin(), u.end(),
    248                                               double(0)) / u.size();
    249         double var = 0;
    250         double skew = 0;
    251         double kurtosis = 0;
    252         for (int i = 0; i < u.size(); ++i)
    253         {
    254             double d = (u[i] - mean);
    255             double d2 = sqr(d);
    256             var += d2;
    257             skew += d * d2;
    258             kurtosis += d2 * d2;
    259         }
    260         var /= u.size();
    261         double dev = std::sqrt(var);
    262         skew /= u.size() * dev * var;
    263         kurtosis /= u.size() * var * var;
    264         kurtosis -= 3;
    265         double x_mean = (1 - d.p()) / d.p();
    266         double x_var = x_mean / d.p();
    267         double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
    268         double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
    269         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
    270         assert(std::abs((var - x_var) / x_var) < 0.01);
    271         assert(std::abs((skew - x_skew) / x_skew) < 0.01);
    272         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
    273     }
    274 }
    275