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 void
     33 test1()
     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 (unsigned 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.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 void
     76 test2()
     77 {
     78     typedef std::geometric_distribution<> D;
     79     typedef std::mt19937 G;
     80     G g;
     81     D d(0.05);
     82     const int N = 1000000;
     83     std::vector<D::result_type> u;
     84     for (int i = 0; i < N; ++i)
     85     {
     86         D::result_type v = d(g);
     87         assert(d.min() <= v && v <= d.max());
     88         u.push_back(v);
     89     }
     90     double mean = std::accumulate(u.begin(), u.end(),
     91                                           double(0)) / u.size();
     92     double var = 0;
     93     double skew = 0;
     94     double kurtosis = 0;
     95     for (unsigned i = 0; i < u.size(); ++i)
     96     {
     97         double dbl = (u[i] - mean);
     98         double d2 = sqr(dbl);
     99         var += d2;
    100         skew += dbl * d2;
    101         kurtosis += d2 * d2;
    102     }
    103     var /= u.size();
    104     double dev = std::sqrt(var);
    105     skew /= u.size() * dev * var;
    106     kurtosis /= u.size() * var * var;
    107     kurtosis -= 3;
    108     double x_mean = (1 - d.p()) / d.p();
    109     double x_var = x_mean / d.p();
    110     double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
    111     double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
    112     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
    113     assert(std::abs((var - x_var) / x_var) < 0.01);
    114     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
    115     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
    116 }
    117 
    118 void
    119 test3()
    120 {
    121     typedef std::geometric_distribution<> D;
    122     typedef std::minstd_rand G;
    123     G g;
    124     D d(.25);
    125     const int N = 1000000;
    126     std::vector<D::result_type> u;
    127     for (int i = 0; i < N; ++i)
    128     {
    129         D::result_type v = d(g);
    130         assert(d.min() <= v && v <= d.max());
    131         u.push_back(v);
    132     }
    133     double mean = std::accumulate(u.begin(), u.end(),
    134                                           double(0)) / u.size();
    135     double var = 0;
    136     double skew = 0;
    137     double kurtosis = 0;
    138     for (unsigned i = 0; i < u.size(); ++i)
    139     {
    140         double dbl = (u[i] - mean);
    141         double d2 = sqr(dbl);
    142         var += d2;
    143         skew += dbl * d2;
    144         kurtosis += d2 * d2;
    145     }
    146     var /= u.size();
    147     double dev = std::sqrt(var);
    148     skew /= u.size() * dev * var;
    149     kurtosis /= u.size() * var * var;
    150     kurtosis -= 3;
    151     double x_mean = (1 - d.p()) / d.p();
    152     double x_var = x_mean / d.p();
    153     double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
    154     double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
    155     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
    156     assert(std::abs((var - x_var) / x_var) < 0.01);
    157     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
    158     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
    159 }
    160 
    161 void
    162 test4()
    163 {
    164     typedef std::geometric_distribution<> D;
    165     typedef std::mt19937 G;
    166     G g;
    167     D d(0.5);
    168     const int N = 1000000;
    169     std::vector<D::result_type> u;
    170     for (int i = 0; i < N; ++i)
    171     {
    172         D::result_type v = d(g);
    173         assert(d.min() <= v && v <= d.max());
    174         u.push_back(v);
    175     }
    176     double mean = std::accumulate(u.begin(), u.end(),
    177                                           double(0)) / u.size();
    178     double var = 0;
    179     double skew = 0;
    180     double kurtosis = 0;
    181     for (unsigned i = 0; i < u.size(); ++i)
    182     {
    183         double dbl = (u[i] - mean);
    184         double d2 = sqr(dbl);
    185         var += d2;
    186         skew += dbl * d2;
    187         kurtosis += d2 * d2;
    188     }
    189     var /= u.size();
    190     double dev = std::sqrt(var);
    191     skew /= u.size() * dev * var;
    192     kurtosis /= u.size() * var * var;
    193     kurtosis -= 3;
    194     double x_mean = (1 - d.p()) / d.p();
    195     double x_var = x_mean / d.p();
    196     double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
    197     double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
    198     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
    199     assert(std::abs((var - x_var) / x_var) < 0.01);
    200     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
    201     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
    202 }
    203 
    204 void
    205 test5()
    206 {
    207     typedef std::geometric_distribution<> D;
    208     typedef std::mt19937 G;
    209     G g;
    210     D d(0.75);
    211     const int N = 1000000;
    212     std::vector<D::result_type> u;
    213     for (int i = 0; i < N; ++i)
    214     {
    215         D::result_type v = d(g);
    216         assert(d.min() <= v && v <= d.max());
    217         u.push_back(v);
    218     }
    219     double mean = std::accumulate(u.begin(), u.end(),
    220                                           double(0)) / u.size();
    221     double var = 0;
    222     double skew = 0;
    223     double kurtosis = 0;
    224     for (unsigned i = 0; i < u.size(); ++i)
    225     {
    226         double dbl = (u[i] - mean);
    227         double d2 = sqr(dbl);
    228         var += d2;
    229         skew += dbl * d2;
    230         kurtosis += d2 * d2;
    231     }
    232     var /= u.size();
    233     double dev = std::sqrt(var);
    234     skew /= u.size() * dev * var;
    235     kurtosis /= u.size() * var * var;
    236     kurtosis -= 3;
    237     double x_mean = (1 - d.p()) / d.p();
    238     double x_var = x_mean / d.p();
    239     double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
    240     double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
    241     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
    242     assert(std::abs((var - x_var) / x_var) < 0.01);
    243     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
    244     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
    245 }
    246 
    247 void
    248 test6()
    249 {
    250     typedef std::geometric_distribution<> D;
    251     typedef std::mt19937 G;
    252     G g;
    253     D d(0.96875);
    254     const int N = 1000000;
    255     std::vector<D::result_type> u;
    256     for (int i = 0; i < N; ++i)
    257     {
    258         D::result_type v = d(g);
    259         assert(d.min() <= v && v <= d.max());
    260         u.push_back(v);
    261     }
    262     double mean = std::accumulate(u.begin(), u.end(),
    263                                           double(0)) / u.size();
    264     double var = 0;
    265     double skew = 0;
    266     double kurtosis = 0;
    267     for (unsigned i = 0; i < u.size(); ++i)
    268     {
    269         double dbl = (u[i] - mean);
    270         double d2 = sqr(dbl);
    271         var += d2;
    272         skew += dbl * d2;
    273         kurtosis += d2 * d2;
    274     }
    275     var /= u.size();
    276     double dev = std::sqrt(var);
    277     skew /= u.size() * dev * var;
    278     kurtosis /= u.size() * var * var;
    279     kurtosis -= 3;
    280     double x_mean = (1 - d.p()) / d.p();
    281     double x_var = x_mean / d.p();
    282     double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
    283     double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
    284     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
    285     assert(std::abs((var - x_var) / x_var) < 0.01);
    286     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
    287     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
    288 }
    289 
    290 int main()
    291 {
    292     test1();
    293     test2();
    294     test3();
    295     test4();
    296     test5();
    297     test6();
    298 }
    299