157 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 exponential_distribution15 16// template<class _URNG> result_type operator()(_URNG& g);17 18#include <random>19#include <cassert>20#include <cmath>21#include <cstddef>22#include <numeric>23#include <vector>24 25#include "test_macros.h"26 27template <class T>28inline29T30sqr(T x)31{32 return x * x;33}34 35int main(int, char**)36{37 {38 typedef std::exponential_distribution<> D;39 typedef std::mt19937 G;40 G g;41 D d(.75);42 const int N = 1000000;43 std::vector<D::result_type> u;44 for (int i = 0; i < N; ++i)45 {46 D::result_type v = d(g);47 assert(d.min() < v);48 u.push_back(v);49 }50 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();51 double var = 0;52 double skew = 0;53 double kurtosis = 0;54 for (std::size_t i = 0; i < u.size(); ++i)55 {56 double dbl = (u[i] - mean);57 double d2 = sqr(dbl);58 var += d2;59 skew += dbl * d2;60 kurtosis += d2 * d2;61 }62 var /= u.size();63 double dev = std::sqrt(var);64 skew /= u.size() * dev * var;65 kurtosis /= u.size() * var * var;66 kurtosis -= 3;67 double x_mean = 1/d.lambda();68 double x_var = 1/sqr(d.lambda());69 double x_skew = 2;70 double x_kurtosis = 6;71 assert(std::abs((mean - x_mean) / x_mean) < 0.01);72 assert(std::abs((var - x_var) / x_var) < 0.01);73 assert(std::abs((skew - x_skew) / x_skew) < 0.01);74 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);75 }76 {77 typedef std::exponential_distribution<> D;78 typedef std::mt19937 G;79 G g;80 D d(1);81 const int N = 1000000;82 std::vector<D::result_type> u;83 for (int i = 0; i < N; ++i)84 {85 D::result_type v = d(g);86 assert(d.min() < v);87 u.push_back(v);88 }89 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();90 double var = 0;91 double skew = 0;92 double kurtosis = 0;93 for (std::size_t i = 0; i < u.size(); ++i)94 {95 double dbl = (u[i] - mean);96 double d2 = sqr(dbl);97 var += d2;98 skew += dbl * d2;99 kurtosis += d2 * d2;100 }101 var /= u.size();102 double dev = std::sqrt(var);103 skew /= u.size() * dev * var;104 kurtosis /= u.size() * var * var;105 kurtosis -= 3;106 double x_mean = 1/d.lambda();107 double x_var = 1/sqr(d.lambda());108 double x_skew = 2;109 double x_kurtosis = 6;110 assert(std::abs((mean - x_mean) / x_mean) < 0.01);111 assert(std::abs((var - x_var) / x_var) < 0.01);112 assert(std::abs((skew - x_skew) / x_skew) < 0.01);113 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);114 }115 {116 typedef std::exponential_distribution<> D;117 typedef std::mt19937 G;118 G g;119 D d(10);120 const int N = 1000000;121 std::vector<D::result_type> u;122 for (int i = 0; i < N; ++i)123 {124 D::result_type v = d(g);125 assert(d.min() < v);126 u.push_back(v);127 }128 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();129 double var = 0;130 double skew = 0;131 double kurtosis = 0;132 for (std::size_t i = 0; i < u.size(); ++i)133 {134 double dbl = (u[i] - mean);135 double d2 = sqr(dbl);136 var += d2;137 skew += dbl * d2;138 kurtosis += d2 * d2;139 }140 var /= u.size();141 double dev = std::sqrt(var);142 skew /= u.size() * dev * var;143 kurtosis /= u.size() * var * var;144 kurtosis -= 3;145 double x_mean = 1/d.lambda();146 double x_var = 1/sqr(d.lambda());147 double x_skew = 2;148 double x_kurtosis = 6;149 assert(std::abs((mean - x_mean) / x_mean) < 0.01);150 assert(std::abs((var - x_var) / x_var) < 0.01);151 assert(std::abs((skew - x_skew) / x_skew) < 0.01);152 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);153 }154 155 return 0;156}157