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