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