test/std/numerics/rand/rand.dis/rand.dist.norm/rand.dist.norm.f/eval.pass.cpp

   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 D::param_type P;
  55         typedef std::mt19937 G;
  56         G g;
  57         D d(2, 4);
  58         const int N = 100000;
  59         std::vector<D::result_type> u;
  60         for (int i = 0; i < N; ++i)
  61         {
  62             D::result_type v = d(g);
  63             assert(v >= 0);
  64             u.push_back(v);
  65         }
  66         std::sort(u.begin(), u.end());
  67         for (int i = 0; i < N; ++i)
  68             assert(std::abs(f(u[i], d.m(), d.n()) - double(i)/N) < .01);
  69     }
  70     {
  71         typedef std::fisher_f_distribution<> D;
  72         typedef D::param_type P;
  73         typedef std::mt19937 G;
  74         G g;
  75         D d(4, 2);
  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);
  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], d.m(), d.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         const int N = 100000;
  95         std::vector<D::result_type> u;
  96         for (int i = 0; i < N; ++i)
  97         {
  98             D::result_type v = d(g);
  99             assert(v >= 0);
 100             u.push_back(v);
 101         }
 102         std::sort(u.begin(), u.end());
 103         for (int i = 0; i < N; ++i)
 104             assert(std::abs(f(u[i], d.m(), d.n()) - double(i)/N) < .01);
 105     }
 106 }