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

source code of libcxx/test/std/numerics/rand/rand.dist/rand.dist.pois/rand.dist.pois.gamma/eval.pass.cpp