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