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 RealType = double>
13 // class piecewise_constant_distribution
14
15 // template<class _URNG> result_type operator()(_URNG& g, const param_type& parm);
16
17 #include <random>
18 #include <vector>
19 #include <iterator>
20 #include <numeric>
21 #include <cassert>
22
23 template <class T>
24 inline
25 T
sqr(T x)26 sqr(T x)
27 {
28 return x*x;
29 }
30
main()31 int main()
32 {
33 {
34 typedef std::piecewise_constant_distribution<> D;
35 typedef D::param_type P;
36 typedef std::mt19937_64 G;
37 G g;
38 double b[] = {10, 14, 16, 17};
39 double p[] = {25, 62.5, 12.5};
40 const size_t Np = sizeof(p) / sizeof(p[0]);
41 D d;
42 P pa(b, b+Np+1, p);
43 const int N = 1000000;
44 std::vector<D::result_type> u;
45 for (int i = 0; i < N; ++i)
46 {
47 D::result_type v = d(g, pa);
48 assert(10 <= v && v < 17);
49 u.push_back(v);
50 }
51 std::vector<double> prob(std::begin(p), std::end(p));
52 double s = std::accumulate(prob.begin(), prob.end(), 0.0);
53 for (int i = 0; i < prob.size(); ++i)
54 prob[i] /= s;
55 std::sort(u.begin(), u.end());
56 for (int i = 0; i < Np; ++i)
57 {
58 typedef std::vector<D::result_type>::iterator I;
59 I lb = std::lower_bound(u.begin(), u.end(), b[i]);
60 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
61 const size_t Ni = ub - lb;
62 if (prob[i] == 0)
63 assert(Ni == 0);
64 else
65 {
66 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
67 double mean = std::accumulate(lb, ub, 0.0) / Ni;
68 double var = 0;
69 double skew = 0;
70 double kurtosis = 0;
71 for (I j = lb; j != ub; ++j)
72 {
73 double d = (*j - mean);
74 double d2 = sqr(d);
75 var += d2;
76 skew += d * d2;
77 kurtosis += d2 * d2;
78 }
79 var /= Ni;
80 double dev = std::sqrt(var);
81 skew /= Ni * dev * var;
82 kurtosis /= Ni * var * var;
83 kurtosis -= 3;
84 double x_mean = (b[i+1] + b[i]) / 2;
85 double x_var = sqr(b[i+1] - b[i]) / 12;
86 double x_skew = 0;
87 double x_kurtosis = -6./5;
88 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
89 assert(std::abs((var - x_var) / x_var) < 0.01);
90 assert(std::abs(skew - x_skew) < 0.01);
91 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
92 }
93 }
94 }
95 }
96