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 weibull_distribution
15
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>
28inline
29T
30sqr(T x)
31{
32 return x * x;
33}
34
35int main(int, char**)
36{
37 {
38 typedef std::weibull_distribution<> D;
39 typedef std::mt19937 G;
40 G g;
41 D d(0.5, 2);
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(x: 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(x: var);
64 skew /= u.size() * dev * var;
65 kurtosis /= u.size() * var * var;
66 kurtosis -= 3;
67 double x_mean = d.b() * std::tgamma(1 + 1/d.a());
68 double x_var = sqr(d.b()) * std::tgamma(1 + 2/d.a()) - sqr(x_mean);
69 double x_skew = (sqr(d.b())*d.b() * std::tgamma(1 + 3/d.a()) -
70 3*x_mean*x_var - sqr(x_mean)*x_mean) /
71 (std::sqrt(x: x_var)*x_var);
72 double x_kurtosis = (sqr(sqr(d.b())) * std::tgamma(1 + 4/d.a()) -
73 4*x_skew*x_var*sqrt(x: x_var)*x_mean -
74 6*sqr(x_mean)*x_var - sqr(sqr(x_mean))) / sqr(x_var) - 3;
75 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
76 assert(std::abs((var - x_var) / x_var) < 0.01);
77 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
78 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
79 }
80 {
81 typedef std::weibull_distribution<> D;
82 typedef std::mt19937 G;
83 G g;
84 D d(1, .5);
85 const int N = 1000000;
86 std::vector<D::result_type> u;
87 for (int i = 0; i < N; ++i)
88 {
89 D::result_type v = d(g);
90 assert(d.min() <= v);
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 (std::size_t 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(x: var);
107 skew /= u.size() * dev * var;
108 kurtosis /= u.size() * var * var;
109 kurtosis -= 3;
110 double x_mean = d.b() * std::tgamma(1 + 1/d.a());
111 double x_var = sqr(d.b()) * std::tgamma(1 + 2/d.a()) - sqr(x_mean);
112 double x_skew = (sqr(d.b())*d.b() * std::tgamma(1 + 3/d.a()) -
113 3*x_mean*x_var - sqr(x_mean)*x_mean) /
114 (std::sqrt(x: x_var)*x_var);
115 double x_kurtosis = (sqr(sqr(d.b())) * std::tgamma(1 + 4/d.a()) -
116 4*x_skew*x_var*sqrt(x: x_var)*x_mean -
117 6*sqr(x_mean)*x_var - sqr(sqr(x_mean))) / sqr(x_var) - 3;
118 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
119 assert(std::abs((var - x_var) / x_var) < 0.01);
120 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
121 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
122 }
123 {
124 typedef std::weibull_distribution<> D;
125 typedef std::mt19937 G;
126 G g;
127 D d(2, 3);
128 const int N = 1000000;
129 std::vector<D::result_type> u;
130 for (int i = 0; i < N; ++i)
131 {
132 D::result_type v = d(g);
133 assert(d.min() <= v);
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 (std::size_t 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(x: var);
150 skew /= u.size() * dev * var;
151 kurtosis /= u.size() * var * var;
152 kurtosis -= 3;
153 double x_mean = d.b() * std::tgamma(1 + 1/d.a());
154 double x_var = sqr(d.b()) * std::tgamma(1 + 2/d.a()) - sqr(x_mean);
155 double x_skew = (sqr(d.b())*d.b() * std::tgamma(1 + 3/d.a()) -
156 3*x_mean*x_var - sqr(x_mean)*x_mean) /
157 (std::sqrt(x: x_var)*x_var);
158 double x_kurtosis = (sqr(sqr(d.b())) * std::tgamma(1 + 4/d.a()) -
159 4*x_skew*x_var*sqrt(x: x_var)*x_mean -
160 6*sqr(x_mean)*x_var - sqr(sqr(x_mean))) / sqr(x_var) - 3;
161 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
162 assert(std::abs((var - x_var) / x_var) < 0.01);
163 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
164 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
165 }
166
167 return 0;
168}
169

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