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 // REQUIRES: long_tests 11 12 // <random> 13 14 // template<class RealType = double> 15 // class fisher_f_distribution 16 17 // template<class _URNG> result_type operator()(_URNG& g); 18 19 #include <random> 20 #include <cassert> 21 #include <vector> 22 #include <algorithm> 23 #include <cmath> 24 25 double fac(double x) 26 { 27 double r = 1; 28 for (; x > 1; --x) 29 r *= x; 30 return r; 31 } 32 33 double 34 I(double x, unsigned a, unsigned b) 35 { 36 double r = 0; 37 for (int j = a; static_cast<unsigned>(j) <= a+b-1; ++j) 38 r += fac(a+b-1)/(fac(j) * fac(a + b - 1 - j)) * std::pow(x, j) * 39 std::pow(1-x, a+b-1-j); 40 return r; 41 } 42 43 double 44 f(double x, double m, double n) 45 { 46 return I(m * x / (m*x + n), static_cast<unsigned>(m/2), static_cast<unsigned>(n/2)); 47 } 48 49 int main() 50 { 51 // Purposefully only testing even integral values of m and n (for now) 52 { 53 typedef std::fisher_f_distribution<> D; 54 typedef std::mt19937 G; 55 G g; 56 D d(2, 4); 57 const int N = 100000; 58 std::vector<D::result_type> u; 59 for (int i = 0; i < N; ++i) 60 { 61 D::result_type v = d(g); 62 assert(v >= 0); 63 u.push_back(v); 64 } 65 std::sort(u.begin(), u.end()); 66 for (int i = 0; i < N; ++i) 67 assert(std::abs(f(u[i], d.m(), d.n()) - double(i)/N) < .01); 68 } 69 { 70 typedef std::fisher_f_distribution<> D; 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 std::mt19937 G; 89 G g; 90 D d(18, 20); 91 const int N = 100000; 92 std::vector<D::result_type> u; 93 for (int i = 0; i < N; ++i) 94 { 95 D::result_type v = d(g); 96 assert(v >= 0); 97 u.push_back(v); 98 } 99 std::sort(u.begin(), u.end()); 100 for (int i = 0; i < N; ++i) 101 assert(std::abs(f(u[i], d.m(), d.n()) - double(i)/N) < .01); 102 } 103 } 104
