113 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// <random>10 11// class bernoulli_distribution12 13// template<class _URNG> result_type operator()(_URNG& g, const param_type& parm);14 15#include <random>16#include <cassert>17#include <cmath>18#include <cstddef>19#include <numeric>20#include <vector>21 22#include "test_macros.h"23 24template <class T>25inline26T27sqr(T x)28{29 return x * x;30}31 32int main(int, char**)33{34 {35 typedef std::bernoulli_distribution D;36 typedef D::param_type P;37 typedef std::minstd_rand G;38 G g;39 D d(.75);40 P p(.25);41 const int N = 100000;42 std::vector<D::result_type> u;43 for (int i = 0; i < N; ++i)44 u.push_back(d(g, p));45 double mean = std::accumulate(u.begin(), u.end(),46 double(0)) / u.size();47 double var = 0;48 double skew = 0;49 double kurtosis = 0;50 for (std::size_t i = 0; i < u.size(); ++i)51 {52 double dbl = (u[i] - mean);53 double d2 = sqr(dbl);54 var += d2;55 skew += dbl * d2;56 kurtosis += d2 * d2;57 }58 var /= u.size();59 double dev = std::sqrt(var);60 skew /= u.size() * dev * var;61 kurtosis /= u.size() * var * var;62 kurtosis -= 3;63 double x_mean = p.p();64 double x_var = p.p()*(1-p.p());65 double x_skew = (1 - 2 * p.p())/std::sqrt(x_var);66 double x_kurtosis = (6 * sqr(p.p()) - 6 * p.p() + 1)/x_var;67 assert(std::abs((mean - x_mean) / x_mean) < 0.01);68 assert(std::abs((var - x_var) / x_var) < 0.01);69 assert(std::abs((skew - x_skew) / x_skew) < 0.02);70 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.05);71 }72 {73 typedef std::bernoulli_distribution D;74 typedef D::param_type P;75 typedef std::minstd_rand G;76 G g;77 D d(.25);78 P p(.75);79 const int N = 100000;80 std::vector<D::result_type> u;81 for (int i = 0; i < N; ++i)82 u.push_back(d(g, p));83 double mean = std::accumulate(u.begin(), u.end(),84 double(0)) / u.size();85 double var = 0;86 double skew = 0;87 double kurtosis = 0;88 for (std::size_t i = 0; i < u.size(); ++i)89 {90 double dbl = (u[i] - mean);91 double d2 = sqr(dbl);92 var += d2;93 skew += dbl * d2;94 kurtosis += d2 * d2;95 }96 var /= u.size();97 double dev = std::sqrt(var);98 skew /= u.size() * dev * var;99 kurtosis /= u.size() * var * var;100 kurtosis -= 3;101 double x_mean = p.p();102 double x_var = p.p()*(1-p.p());103 double x_skew = (1 - 2 * p.p())/std::sqrt(x_var);104 double x_kurtosis = (6 * sqr(p.p()) - 6 * p.p() + 1)/x_var;105 assert(std::abs((mean - x_mean) / x_mean) < 0.01);106 assert(std::abs((var - x_var) / x_var) < 0.01);107 assert(std::abs((skew - x_skew) / x_skew) < 0.02);108 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.05);109 }110 111 return 0;112}113