////////////////////////////////////////////////////////////////////////////////
// File: f_distribution.c                                                     //
// Routine(s):                                                                //
//    F_Distribution                                                          //
////////////////////////////////////////////////////////////////////////////////
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
//                         Externally Defined Routines                        //
extern double Beta_Distribution(double x, double a, double b);

////////////////////////////////////////////////////////////////////////////////
// double F_Distribution( double x, int v1, int v2 )                          //
//                                                                            //
//  Description:                                                              //
//     If X1 and X2 are independent chi-square distributed random variables   //
//     with v1 and v2 degrees of freedom respectively, then the random        //
//     variable F = (X1/v1) / (X2/v2) has an F-distribution with v1 and v2    //
//     degrees of freedom.                                                    //
//                                                                            //
//     The F-distribution is the Pr[F < x] which equals the integral from     //
//     -inf to x of the density                                               //
//                               0                           if f < 0,        //
//       [v1^(v1/2) v2^(v2/2) / B(v1/2,v2/2)]                                 //
//                 * f^(v1/2-1) * (v2+v1 f)^(-(v1+v2)/2))   if f >= 0,        //
//     where v1 >= 1, v2 >= 1, and B(,) is the (complete) beta function.      //
//                                                                            //
//     By making the change of variables: g = v1*f / (v2 + v1*f),             //
//                   F(x,v1,v2) = B(v1*x / (v2 + v1*x), v1/2, v2/2),          //
//     where B(,,) is the incomplete beta function.                           //
//                                                                            //
//  Arguments:                                                                //
//     double x   The upper limit of the integral of the density given above. //
//     int   v1   The number of degrees of freedom of the numerator of the    //
//                F-test, i.e. (v1/2 - 1) is the exponent of f in the         //
//                integrand above.  Note v1 >= 1.                             //
//     int   v2   The number of degrees of freedom of the denominator of the  //
//                F-test, i.e. (-(v1+v2)/2) is the exponent of the term       //
//                (v2 + v1 f) in the integrand above.  Note v2 >= 1.          //
//                                                                            //
//  Return Values:                                                            //
//     A real number between 0 and 1.                                         //
//                                                                            //
//  Example:                                                                  //
//     double p, x;                                                           //
//     int    v1, v2;                                                         //
//                                                                            //
//     p = F_Distribution(x, v1, v2);                                         //
////////////////////////////////////////////////////////////////////////////////

double F_Distribution(double f, int v1, int v2)
{
   double a = (double) v1 / 2.0;
   double b = (double) v2 / 2.0;
   double g = a*f;

   if ( f <= 0.0 ) return 0.0;

   return Beta_Distribution( g / (b + g), a, b);
}

double comp_F_Distribution(double f, int v1, int v2)
{
   double a = (double) v1 / 2.0;
   double b = (double) v2 / 2.0;
   double g = a*f;

   if ( f <= 0.0 ) return 1.0;

   return Beta_Distribution( 1.0-g / (b + g), b, a);
}

// Fcache routines to cache the values for -log10(Fprob(1,df)) 
// Fcache_max is 6, which is roughly the upper 1% point, so 99% of pvalues will be retieved from the cache

static double *Fcache=NULL;
static double Fcache_max=6;
static int fc=0;
void init_Fcache( int df ) {

  int i;
  int imax = (int)Fcache_max*10000;
  if ( Fcache != NULL )
    free(Fcache);
  Fcache = (double*)calloc(imax,sizeof(double));

  for(i=0;i<imax;i++)
    Fcache[i] = -log10(comp_F_Distribution(i/10000.0,1, df));
}

double Log10Fpval( double F, double df1, double df2 ) {
  //  if ( !(fc%10000) ) printf("fc %d F %e df1 %e\n", fc, F, df1);
  //  fc++;
  if ( df1 == 1.0 && F < Fcache_max ) {
    //    fc++;
    return Fcache[(int)(10000*F)];
  }
  else {
    return (-log10(comp_F_Distribution(F,(int)df1, (int)df2)));
  }
}