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 fisher_f_distribution
15
16 // template<class _URNG> result_type operator()(_URNG& g, const param_type& parm);
17
18 #include <random>
19 #include <cassert>
20 #include <vector>
21 #include <algorithm>
22 #include <cmath>
23
24 #include "test_macros.h"
25
fac(double x)26 double fac(double x)
27 {
28 double r = 1;
29 for (; x > 1; --x)
30 r *= x;
31 return r;
32 }
33
34 double
I(double x,unsigned a,unsigned b)35 I(double x, unsigned a, unsigned b)
36 {
37 double r = 0;
38 for (int j = a; static_cast<unsigned>(j) <= a+b-1; ++j)
39 r += fac(a+b-1)/(fac(j) * fac(a + b - 1 - j)) * std::pow(x, j) *
40 std::pow(1-x, a+b-1-j);
41 return r;
42 }
43
44 double
f(double x,double m,double n)45 f(double x, double m, double n)
46 {
47 return I(m * x / (m*x + n), static_cast<unsigned>(m/2), static_cast<unsigned>(n/2));
48 }
49
main(int,char **)50 int main(int, char**)
51 {
52 // Purposefully only testing even integral values of m and n (for now)
53 {
54 typedef std::fisher_f_distribution<> D;
55 typedef D::param_type P;
56 typedef std::mt19937 G;
57 G g;
58 D d(2, 4);
59 P p(4, 2);
60 const int N = 100000;
61 std::vector<D::result_type> u;
62 for (int i = 0; i < N; ++i)
63 {
64 D::result_type v = d(g, p);
65 assert(v >= 0);
66 u.push_back(v);
67 }
68 std::sort(u.begin(), u.end());
69 for (int i = 0; i < N; ++i)
70 assert(std::abs(f(u[i], p.m(), p.n()) - double(i)/N) < .01);
71 }
72 {
73 typedef std::fisher_f_distribution<> D;
74 typedef D::param_type P;
75 typedef std::mt19937 G;
76 G g;
77 D d(4, 2);
78 P p(6, 8);
79 const int N = 100000;
80 std::vector<D::result_type> u;
81 for (int i = 0; i < N; ++i)
82 {
83 D::result_type v = d(g, p);
84 assert(v >= 0);
85 u.push_back(v);
86 }
87 std::sort(u.begin(), u.end());
88 for (int i = 0; i < N; ++i)
89 assert(std::abs(f(u[i], p.m(), p.n()) - double(i)/N) < .01);
90 }
91 {
92 typedef std::fisher_f_distribution<> D;
93 typedef D::param_type P;
94 typedef std::mt19937 G;
95 G g;
96 D d(18, 20);
97 P p(16, 14);
98 const int N = 100000;
99 std::vector<D::result_type> u;
100 for (int i = 0; i < N; ++i)
101 {
102 D::result_type v = d(g, p);
103 assert(v >= 0);
104 u.push_back(v);
105 }
106 std::sort(u.begin(), u.end());
107 for (int i = 0; i < N; ++i)
108 assert(std::abs(f(u[i], p.m(), p.n()) - double(i)/N) < .01);
109 }
110
111 return 0;
112 }
113