Home | History | Annotate | Download | only in rand.dist.norm.f 
      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, const param_type& parm);
     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         P p(4, 2);
     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, p);
     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], p.m(), p.n()) - double(i)/N) < .01);
     68     }
     69     {
     70         typedef std::fisher_f_distribution<> D;
     71         typedef D::param_type P;
     72         typedef std::mt19937 G;
     73         G g;
     74         D d(4, 2);
     75         P p(6, 8);
     76         const int N = 100000;
     77         std::vector<D::result_type> u;
     78         for (int i = 0; i < N; ++i)
     79         {
     80             D::result_type v = d(g, p);
     81             assert(v >= 0);
     82             u.push_back(v);
     83         }
     84         std::sort(u.begin(), u.end());
     85         for (int i = 0; i < N; ++i)
     86             assert(std::abs(f(u[i], p.m(), p.n()) - double(i)/N) < .01);
     87     }
     88     {
     89         typedef std::fisher_f_distribution<> D;
     90         typedef D::param_type P;
     91         typedef std::mt19937 G;
     92         G g;
     93         D d(18, 20);
     94         P p(16, 14);
     95         const int N = 100000;
     96         std::vector<D::result_type> u;
     97         for (int i = 0; i < N; ++i)
     98         {
     99             D::result_type v = d(g, p);
    100             assert(v >= 0);
    101             u.push_back(v);
    102         }
    103         std::sort(u.begin(), u.end());
    104         for (int i = 0; i < N; ++i)
    105             assert(std::abs(f(u[i], p.m(), p.n()) - double(i)/N) < .01);
    106     }
    107 }
    108