• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 poisson_distribution
16 
17 // template<class _URNG> result_type operator()(_URNG& g, const param_type& parm);
18 
19 #include <random>
20 #include <cassert>
21 #include <vector>
22 #include <numeric>
23 
24 template <class T>
25 inline
26 T
sqr(T x)27 sqr(T x)
28 {
29     return x * x;
30 }
31 
main()32 int main()
33 {
34     {
35         typedef std::poisson_distribution<> D;
36         typedef D::param_type P;
37         typedef std::minstd_rand G;
38         G g;
39         D d(.75);
40         P p(2);
41         const int N = 100000;
42         std::vector<double> u;
43         for (int i = 0; i < N; ++i)
44         {
45             D::result_type v = d(g, p);
46             assert(d.min() <= v && v <= d.max());
47             u.push_back(v);
48         }
49         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
50         double var = 0;
51         double skew = 0;
52         double kurtosis = 0;
53         for (int i = 0; i < u.size(); ++i)
54         {
55             double d = (u[i] - mean);
56             double d2 = sqr(d);
57             var += d2;
58             skew += d * d2;
59             kurtosis += d2 * d2;
60         }
61         var /= u.size();
62         double dev = std::sqrt(var);
63         skew /= u.size() * dev * var;
64         kurtosis /= u.size() * var * var;
65         kurtosis -= 3;
66         double x_mean = p.mean();
67         double x_var = p.mean();
68         double x_skew = 1 / std::sqrt(x_var);
69         double x_kurtosis = 1 / x_var;
70         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
71         assert(std::abs((var - x_var) / x_var) < 0.01);
72         assert(std::abs((skew - x_skew) / x_skew) < 0.01);
73         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
74     }
75     {
76         typedef std::poisson_distribution<> D;
77         typedef D::param_type P;
78         typedef std::minstd_rand G;
79         G g;
80         D d(2);
81         P p(.75);
82         const int N = 100000;
83         std::vector<double> u;
84         for (int i = 0; i < N; ++i)
85         {
86             D::result_type v = d(g, p);
87             assert(d.min() <= v && v <= d.max());
88             u.push_back(v);
89         }
90         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
91         double var = 0;
92         double skew = 0;
93         double kurtosis = 0;
94         for (int i = 0; i < u.size(); ++i)
95         {
96             double d = (u[i] - mean);
97             double d2 = sqr(d);
98             var += d2;
99             skew += d * d2;
100             kurtosis += d2 * d2;
101         }
102         var /= u.size();
103         double dev = std::sqrt(var);
104         skew /= u.size() * dev * var;
105         kurtosis /= u.size() * var * var;
106         kurtosis -= 3;
107         double x_mean = p.mean();
108         double x_var = p.mean();
109         double x_skew = 1 / std::sqrt(x_var);
110         double x_kurtosis = 1 / x_var;
111         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
112         assert(std::abs((var - x_var) / x_var) < 0.01);
113         assert(std::abs((skew - x_skew) / x_skew) < 0.01);
114         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.04);
115     }
116     {
117         typedef std::poisson_distribution<> D;
118         typedef D::param_type P;
119         typedef std::mt19937 G;
120         G g;
121         D d(2);
122         P p(20);
123         const int N = 1000000;
124         std::vector<double> u;
125         for (int i = 0; i < N; ++i)
126         {
127             D::result_type v = d(g, p);
128             assert(d.min() <= v && v <= d.max());
129             u.push_back(v);
130         }
131         double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
132         double var = 0;
133         double skew = 0;
134         double kurtosis = 0;
135         for (int i = 0; i < u.size(); ++i)
136         {
137             double d = (u[i] - mean);
138             double d2 = sqr(d);
139             var += d2;
140             skew += d * d2;
141             kurtosis += d2 * d2;
142         }
143         var /= u.size();
144         double dev = std::sqrt(var);
145         skew /= u.size() * dev * var;
146         kurtosis /= u.size() * var * var;
147         kurtosis -= 3;
148         double x_mean = p.mean();
149         double x_var = p.mean();
150         double x_skew = 1 / std::sqrt(x_var);
151         double x_kurtosis = 1 / x_var;
152         assert(std::abs((mean - x_mean) / x_mean) < 0.01);
153         assert(std::abs((var - x_var) / x_var) < 0.01);
154         assert(std::abs((skew - x_skew) / x_skew) < 0.01);
155         assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
156     }
157 }
158