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