215 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 extreme_value_distribution15 16// template<class _URNG> result_type operator()(_URNG& g, const param_type& parm);17 18#include <random>19#include <cassert>20#include <cmath>21#include <numeric>22#include <vector>23 24#include "test_macros.h"25 26template <class T>27inline28T29sqr(T x)30{31 return x * x;32}33 34void35test1()36{37 typedef std::extreme_value_distribution<> D;38 typedef D::param_type P;39 typedef std::mt19937 G;40 G g;41 D d(-0.5, 1);42 P p(0.5, 2);43 const int N = 1000000;44 std::vector<D::result_type> u;45 for (int i = 0; i < N; ++i)46 {47 D::result_type v = d(g, p);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 (unsigned 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 = p.a() + p.b() * 0.577215665;68 double x_var = sqr(p.b()) * 1.644934067;69 double x_skew = 1.139547;70 double x_kurtosis = 12./5;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 77void78test2()79{80 typedef std::extreme_value_distribution<> D;81 typedef D::param_type P;82 typedef std::mt19937 G;83 G g;84 D d(-0.5, 1);85 P p(1, 2);86 const int N = 1000000;87 std::vector<D::result_type> u;88 for (int i = 0; i < N; ++i)89 {90 D::result_type v = d(g, p);91 u.push_back(v);92 }93 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();94 double var = 0;95 double skew = 0;96 double kurtosis = 0;97 for (unsigned i = 0; i < u.size(); ++i)98 {99 double dbl = (u[i] - mean);100 double d2 = sqr(dbl);101 var += d2;102 skew += dbl * d2;103 kurtosis += d2 * d2;104 }105 var /= u.size();106 double dev = std::sqrt(var);107 skew /= u.size() * dev * var;108 kurtosis /= u.size() * var * var;109 kurtosis -= 3;110 double x_mean = p.a() + p.b() * 0.577215665;111 double x_var = sqr(p.b()) * 1.644934067;112 double x_skew = 1.139547;113 double x_kurtosis = 12./5;114 assert(std::abs((mean - x_mean) / x_mean) < 0.01);115 assert(std::abs((var - x_var) / x_var) < 0.01);116 assert(std::abs((skew - x_skew) / x_skew) < 0.01);117 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);118}119 120void121test3()122{123 typedef std::extreme_value_distribution<> D;124 typedef D::param_type P;125 typedef std::mt19937 G;126 G g;127 D d(-0.5, 1);128 P p(1.5, 3);129 const int N = 1000000;130 std::vector<D::result_type> u;131 for (int i = 0; i < N; ++i)132 {133 D::result_type v = d(g, p);134 u.push_back(v);135 }136 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();137 double var = 0;138 double skew = 0;139 double kurtosis = 0;140 for (unsigned i = 0; i < u.size(); ++i)141 {142 double dbl = (u[i] - mean);143 double d2 = sqr(dbl);144 var += d2;145 skew += dbl * d2;146 kurtosis += d2 * d2;147 }148 var /= u.size();149 double dev = std::sqrt(var);150 skew /= u.size() * dev * var;151 kurtosis /= u.size() * var * var;152 kurtosis -= 3;153 double x_mean = p.a() + p.b() * 0.577215665;154 double x_var = sqr(p.b()) * 1.644934067;155 double x_skew = 1.139547;156 double x_kurtosis = 12./5;157 assert(std::abs((mean - x_mean) / x_mean) < 0.01);158 assert(std::abs((var - x_var) / x_var) < 0.01);159 assert(std::abs((skew - x_skew) / x_skew) < 0.01);160 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);161}162 163void164test4()165{166 typedef std::extreme_value_distribution<> D;167 typedef D::param_type P;168 typedef std::mt19937 G;169 G g;170 D d(-0.5, 1);171 P p(3, 4);172 const int N = 1000000;173 std::vector<D::result_type> u;174 for (int i = 0; i < N; ++i)175 {176 D::result_type v = d(g, p);177 u.push_back(v);178 }179 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();180 double var = 0;181 double skew = 0;182 double kurtosis = 0;183 for (unsigned i = 0; i < u.size(); ++i)184 {185 double dbl = (u[i] - mean);186 double d2 = sqr(dbl);187 var += d2;188 skew += dbl * d2;189 kurtosis += d2 * d2;190 }191 var /= u.size();192 double dev = std::sqrt(var);193 skew /= u.size() * dev * var;194 kurtosis /= u.size() * var * var;195 kurtosis -= 3;196 double x_mean = p.a() + p.b() * 0.577215665;197 double x_var = sqr(p.b()) * 1.644934067;198 double x_skew = 1.139547;199 double x_kurtosis = 12./5;200 assert(std::abs((mean - x_mean) / x_mean) < 0.01);201 assert(std::abs((var - x_var) / x_var) < 0.01);202 assert(std::abs((skew - x_skew) / x_skew) < 0.01);203 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);204}205 206int main(int, char**)207{208 test1();209 test2();210 test3();211 test4();212 213 return 0;214}215