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