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

   1 /* cdf/fdistinv.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 <gsl/gsl_cdf.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25
  26 #include "error.h"
  27
  28 double
  29 gsl_cdf_fdist_Pinv (const double P, const double nu1, const double nu2)
  30 {
  31   double result;
  32   double y;
  33
  34   if (P < 0.0)
  35     {
  36       CDF_ERROR ("P < 0.0", GSL_EDOM);
  37     }
  38   if (P > 1.0)
  39     {
  40       CDF_ERROR ("P > 1.0", GSL_EDOM);
  41     }
  42   if (nu1 < 1.0)
  43     {
  44       CDF_ERROR ("nu1 < 1", GSL_EDOM);
  45     }
  46   if (nu2 < 1.0)
  47     {
  48       CDF_ERROR ("nu2 < 1", GSL_EDOM);
  49     }
  50
  51   if (P < 0.5)
  52     {
  53       y = gsl_cdf_beta_Pinv (P, nu1 / 2.0, nu2 / 2.0);
  54
  55       result = nu2 * y / (nu1 * (1.0 - y));
  56     }
  57   else
  58     {
  59       y = gsl_cdf_beta_Qinv (P, nu2 / 2.0, nu1 / 2.0);
  60
  61       result = nu2 * (1 - y) / (nu1 * y);
  62     }
  63
  64   return result;
  65 }
  66
  67 double
  68 gsl_cdf_fdist_Qinv (const double Q, const double nu1, const double nu2)
  69 {
  70   double result;
  71   double y;
  72
  73   if (Q < 0.0)
  74     {
  75       CDF_ERROR ("Q < 0.0", GSL_EDOM);
  76     }
  77   if (Q > 1.0)
  78     {
  79       CDF_ERROR ("Q > 1.0", GSL_EDOM);
  80     }
  81   if (nu1 < 1.0)
  82     {
  83       CDF_ERROR ("nu1 < 1", GSL_EDOM);
  84     }
  85   if (nu2 < 1.0)
  86     {
  87       CDF_ERROR ("nu2 < 1", GSL_EDOM);
  88     }
  89
  90   if (Q > 0.5)
  91     {
  92       y = gsl_cdf_beta_Qinv (Q, nu1 / 2.0, nu2 / 2.0);
  93
  94       result = nu2 * y / (nu1 * (1.0 - y));
  95     }
  96   else
  97     {
  98       y = gsl_cdf_beta_Pinv (Q, nu2 / 2.0, nu1 / 2.0);
  99
 100       result = nu2 * (1 - y) / (nu1 * y);
 101     }
 102
 103   return result;
 104 }