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