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