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

   1 /* cdf/hypergeometric.c
   2  *
   3  * Copyright (C) 2004 Jason H. Stover.
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
  18  */
  19
  20 /*
  21  * Computes the cumulative distribution function for a hypergeometric
  22  * random variable. A hypergeometric random variable X is the number
  23  * of elements of type 1 in a sample of size t, drawn from a population
  24  * of size n1 + n2, in which n1 are of type 1 and n2 are of type 2.
  25  *
  26  * This algorithm computes Pr( X <= k ) by summing the terms from
  27  * the mass function, Pr( X = k ).
  28  *
  29  * References:
  30  *
  31  * T. Wu. An accurate computation of the hypergeometric distribution
  32  * function. ACM Transactions on Mathematical Software. Volume 19, number 1,
  33  * March 1993.
  34  *  This algorithm is not used, since it requires factoring the
  35  *  numerator and denominator, then cancelling. It is more accurate
  36  *  than the algorithm used here, but the cancellation requires more
  37  *  time than the algorithm used here.
  38  *
  39  * W. Feller. An Introduction to Probability Theory and Its Applications,
  40  * third edition. 1968. Chapter 2, section 6.
  41  */
  42
  43 #include <config.h>
  44 #include <math.h>
  45 #include <gsl/gsl_math.h>
  46 #include <gsl/gsl_errno.h>
  47 #include <gsl/gsl_cdf.h>
  48 #include <gsl/gsl_randist.h>
  49
  50 #include "error.h"
  51
  52 static double
  53 lower_tail (const unsigned int k, const unsigned int n1,
  54             const unsigned int n2, const unsigned int t)
  55 {
  56   double relerr;
  57   int i = k;
  58   double s, P;
  59
  60   s = gsl_ran_hypergeometric_pdf (i, n1, n2, t);
  61   P = s;
  62
  63   while (i > 0)
  64     {
  65       double factor =
  66         (i / (n1 - i + 1.0)) * ((n2 + i - t) / (t - i + 1.0));
  67       s *= factor;
  68       P += s;
  69       relerr = s / P;
  70       if (relerr < GSL_DBL_EPSILON)
  71         break;
  72       i--;
  73     }
  74
  75   return P;
  76 }
  77
  78 static double
  79 upper_tail (const unsigned int k, const unsigned int n1,
  80             const unsigned int n2, const unsigned int t)
  81 {
  82   double relerr;
  83   unsigned int i = k + 1;
  84   double s, Q;
  85
  86   s = gsl_ran_hypergeometric_pdf (i, n1, n2, t);
  87   Q = s;
  88
  89   while (i < t)
  90     {
  91       double factor =
  92         ((n1 - i) / (i + 1.0)) * ((t - i) / (n2 + i + 1.0 - t));
  93       s *= factor;
  94       Q += s;
  95       relerr = s / Q;
  96       if (relerr < GSL_DBL_EPSILON)
  97         break;
  98       i++;
  99     }
 100
 101   return Q;
 102 }
 103
 104
 105
 106
 107 /*
 108  * Pr (X <= k)
 109  */
 110 double
 111 gsl_cdf_hypergeometric_P (const unsigned int k,
 112                           const unsigned int n1,
 113                           const unsigned int n2, const unsigned int t)
 114 {
 115   double P;
 116
 117   if (t > (n1 + n2))
 118     {
 119       CDF_ERROR ("t larger than population size", GSL_EDOM);
 120     }
 121   else if (k >= n1 || k >= t)
 122     {
 123       P = 1.0;
 124     }
 125   else if (k < 0.0)
 126     {
 127       P = 0.0;
 128     }
 129   else
 130     {
 131       double midpoint = ((double)t * n1) / ((double)n1 + (double)n2);
 132
 133       if (k >= midpoint)
 134         {
 135           P = 1 - upper_tail (k, n1, n2, t);
 136         }
 137       else
 138         {
 139           P = lower_tail (k, n1, n2, t);
 140         }
 141     }
 142
 143   return P;
 144 }
 145
 146 /*
 147  * Pr (X > k)
 148  */
 149 double
 150 gsl_cdf_hypergeometric_Q (const unsigned int k,
 151                           const unsigned int n1,
 152                           const unsigned int n2, const unsigned int t)
 153 {
 154   double Q;
 155
 156   if (t > (n1 + n2))
 157     {
 158       CDF_ERROR ("t larger than population size", GSL_EDOM);
 159     }
 160   else if (k >= n1 || k >= t)
 161     {
 162       Q = 0.0;
 163     }
 164   else if (k < 0.0)
 165     {
 166       Q = 1.0;
 167     }
 168   else
 169     {
 170       double midpoint = ((double)t * n1) / ((double)n1 + (double)n2);
 171
 172       if (k < midpoint)
 173         {
 174           Q = 1 - lower_tail (k, n1, n2, t);
 175         }
 176       else
 177         {
 178           Q = upper_tail (k, n1, n2, t);
 179         }
 180     }
 181
 182   return Q;
 183 }