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