1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
|
//===----------------------------------------------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
//
// REQUIRES: long_tests
// <random>
// template<class RealType = double>
// class fisher_f_distribution
// template<class _URNG> result_type operator()(_URNG& g);
#include <random>
#include <cassert>
#include <vector>
#include <algorithm>
#include <cmath>
#include "test_macros.h"
double fac(double x)
{
double r = 1;
for (; x > 1; --x)
r *= x;
return r;
}
double
I(double x, unsigned a, unsigned b)
{
double r = 0;
for (int j = a; static_cast<unsigned>(j) <= a+b-1; ++j)
r += fac(a+b-1)/(fac(j) * fac(a + b - 1 - j)) * std::pow(x, j) *
std::pow(1-x, a+b-1-j);
return r;
}
double
f(double x, double m, double n)
{
return I(m * x / (m*x + n), static_cast<unsigned>(m/2), static_cast<unsigned>(n/2));
}
int main(int, char**)
{
// Purposefully only testing even integral values of m and n (for now)
{
typedef std::fisher_f_distribution<> D;
typedef std::mt19937 G;
G g;
D d(2, 4);
const int N = 100000;
std::vector<D::result_type> u;
for (int i = 0; i < N; ++i)
{
D::result_type v = d(g);
assert(v >= 0);
u.push_back(v);
}
std::sort(u.begin(), u.end());
for (int i = 0; i < N; ++i)
assert(std::abs(f(u[i], d.m(), d.n()) - double(i)/N) < .01);
}
{
typedef std::fisher_f_distribution<> D;
typedef std::mt19937 G;
G g;
D d(4, 2);
const int N = 100000;
std::vector<D::result_type> u;
for (int i = 0; i < N; ++i)
{
D::result_type v = d(g);
assert(v >= 0);
u.push_back(v);
}
std::sort(u.begin(), u.end());
for (int i = 0; i < N; ++i)
assert(std::abs(f(u[i], d.m(), d.n()) - double(i)/N) < .01);
}
{
typedef std::fisher_f_distribution<> D;
typedef std::mt19937 G;
G g;
D d(18, 20);
const int N = 100000;
std::vector<D::result_type> u;
for (int i = 0; i < N; ++i)
{
D::result_type v = d(g);
assert(v >= 0);
u.push_back(v);
}
std::sort(u.begin(), u.end());
for (int i = 0; i < N; ++i)
assert(std::abs(f(u[i], d.m(), d.n()) - double(i)/N) < .01);
}
return 0;
}
|
