1 //===----------------------------------------------------------------------===//
2 //
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4 // See https://llvm.org/LICENSE.txt for license information.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6 //
7 //===----------------------------------------------------------------------===//
8 //
9 // REQUIRES: long_tests
10
11 // <random>
12
13 // template<class RealType = double>
14 // class piecewise_constant_distribution
15
16 // template<class _URNG> result_type operator()(_URNG& g, const param_type& parm);
17
18 #include <random>
19 #include <algorithm>
20 #include <vector>
21 #include <iterator>
22 #include <numeric>
23 #include <cassert>
24 #include <cstddef>
25
26 #include "test_macros.h"
27
28 template <class T>
29 inline
30 T
sqr(T x)31 sqr(T x)
32 {
33 return x*x;
34 }
35
main(int,char **)36 int main(int, char**)
37 {
38 {
39 typedef std::piecewise_constant_distribution<> D;
40 typedef D::param_type P;
41 typedef std::mt19937_64 G;
42 G g;
43 double b[] = {10, 14, 16, 17};
44 double p[] = {25, 62.5, 12.5};
45 const size_t Np = sizeof(p) / sizeof(p[0]);
46 D d;
47 P pa(b, b+Np+1, p);
48 const int N = 1000000;
49 std::vector<D::result_type> u;
50 for (int i = 0; i < N; ++i)
51 {
52 D::result_type v = d(g, pa);
53 assert(10 <= v && v < 17);
54 u.push_back(v);
55 }
56 std::vector<double> prob(std::begin(p), std::end(p));
57 double s = std::accumulate(prob.begin(), prob.end(), 0.0);
58 for (std::size_t i = 0; i < prob.size(); ++i)
59 prob[i] /= s;
60 std::sort(u.begin(), u.end());
61 for (std::size_t i = 0; i < Np; ++i)
62 {
63 typedef std::vector<D::result_type>::iterator I;
64 I lb = std::lower_bound(u.begin(), u.end(), b[i]);
65 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
66 const size_t Ni = ub - lb;
67 if (prob[i] == 0)
68 assert(Ni == 0);
69 else
70 {
71 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
72 double mean = std::accumulate(lb, ub, 0.0) / Ni;
73 double var = 0;
74 double skew = 0;
75 double kurtosis = 0;
76 for (I j = lb; j != ub; ++j)
77 {
78 double dbl = (*j - mean);
79 double d2 = sqr(dbl);
80 var += d2;
81 skew += dbl * d2;
82 kurtosis += d2 * d2;
83 }
84 var /= Ni;
85 double dev = std::sqrt(var);
86 skew /= Ni * dev * var;
87 kurtosis /= Ni * var * var;
88 kurtosis -= 3;
89 double x_mean = (b[i+1] + b[i]) / 2;
90 double x_var = sqr(b[i+1] - b[i]) / 12;
91 double x_skew = 0;
92 double x_kurtosis = -6./5;
93 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
94 assert(std::abs((var - x_var) / x_var) < 0.01);
95 assert(std::abs(skew - x_skew) < 0.01);
96 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
97 }
98 }
99 }
100
101 return 0;
102 }
103