htswanalysis/MACS/lib/gsl/gsl-1.11/cdf/fdist.c

   1 /* cdf/fdist.c
   2  *
   3  * Copyright (C) 2002 Przemyslaw Sliwa and Jason H. Stover.
   4  * Copyright (C) 2006, 2007 Brian Gough
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
  19  */
  20
  21 #include <config.h>
  22 #include <math.h>
  23 #include <gsl/gsl_cdf.h>
  24 #include <gsl/gsl_sf_gamma.h>
  25 #include <gsl/gsl_math.h>
  26 #include <gsl/gsl_errno.h>
  27
  28 #include "error.h"
  29 #include "beta_inc.c"
  30
  31 /*
  32  * Lower tail.
  33  */
  34
  35 double
  36 gsl_cdf_fdist_P (const double x, const double nu1, const double nu2)
  37 {
  38   double P;
  39   double r = nu2 / nu1;
  40
  41   if (x < r)
  42     {
  43       double u = x / (r + x);
  44
  45       P = beta_inc_AXPY (1.0, 0.0, nu1 / 2.0, nu2 / 2.0, u);
  46     }
  47   else
  48     {
  49       double u = r / (r + x);
  50
  51       P = beta_inc_AXPY (-1.0, 1.0, nu2 / 2.0, nu1 / 2.0, u);
  52     }
  53
  54   return P;
  55 }
  56
  57 /*
  58  * Upper tail.
  59  */
  60
  61 double
  62 gsl_cdf_fdist_Q (const double x, const double nu1, const double nu2)
  63 {
  64   double Q;
  65   double r = nu2 / nu1;
  66
  67   if (x < r)
  68     {
  69       double u = x / (r + x);
  70
  71       Q = beta_inc_AXPY (-1.0, 1.0, nu1 / 2.0, nu2 / 2.0, u);
  72     }
  73   else
  74     {
  75       double u = r / (r + x);
  76
  77       Q = beta_inc_AXPY (1.0, 0.0, nu2 / 2.0, nu1 / 2.0, u);
  78     }
  79
  80   return Q;
  81 }