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
25template <class T>
26inline
27T
28sqr(T x)
29{
30    return x * x;
31}
32
33int 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