lib/gsl-extras/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 2 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 0 in a sample of size t, drawn from a population
  24  * of size n1 + n0, in which n1 are of type 1 and n0 are of type 0.
  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 #include <config.h>
  43 #include <math.h>
  44 #include <gsl/gsl_math.h>
  45 #include <gsl/gsl_errno.h>
  46 #include <gsl/gsl_cdf.h>
  47 #include <gsl/gsl_randist.h>
  48 #include "gsl-extras.h"
  49
  50 /*
  51  * Pr (X <= k)
  52  */
  53 double
  54 gslextras_cdf_hypergeometric_P (const unsigned int k,
  55                                 const unsigned int n0,
  56                                 const unsigned int n1,
  57                                 const unsigned int t)
  58 {
  59   unsigned int i;
  60   unsigned int mode;
  61   double P;
  62   double tmp;
  63   double relerr;
  64
  65   if( t > (n0+n1))
  66     {
  67       GSLEXTRAS_CDF_ERROR("t larger than population size",GSL_EDOM);
  68     }
  69   else if( k >= n0 || k >= t)
  70     {
  71       P = 1.0;
  72     }
  73   else if (k < 0.0)
  74     {
  75       P = 0.0;
  76     }
  77   else
  78     {
  79       P = 0.0;
  80       mode = (int) t*n0 / (n0+n1);
  81       relerr = 1.0;
  82       if( k < mode )
  83         {
  84           i = k;
  85           relerr = 1.0;
  86           while(i != UINT_MAX && relerr > GSL_DBL_EPSILON && P < 1.0)
  87             {
  88               tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
  89               P += tmp;
  90               relerr = tmp / P;
  91               i--;
  92             }
  93         }
  94       else
  95         {
  96           i = mode;
  97           relerr = 1.0;
  98           while(i <= k && relerr > GSL_DBL_EPSILON && P < 1.0)
  99             {
 100               tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
 101               P += tmp;
 102               relerr = tmp / P;
 103               i++;
 104             }
 105           i = mode - 1;
 106           relerr = 1.0;
 107           while( i != UINT_MAX && relerr > GSL_DBL_EPSILON && P < 1.0)
 108             {
 109               tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
 110               P += tmp;
 111               relerr = tmp / P;
 112               i--;
 113             }
 114         }
 115       /*
 116        * Hack to get rid of a pesky error when the sum
 117        * gets slightly above 1.0.
 118        */
 119       P = GSL_MIN_DBL (P, 1.0);
 120     }
 121   return P;
 122 }
 123
 124 /*
 125  * Pr (X > k)
 126  */
 127 double
 128 gslextras_cdf_hypergeometric_Q (const unsigned int k,
 129                                 const unsigned int n0,
 130                                 const unsigned int n1,
 131                                 const unsigned int t)
 132 {
 133   unsigned int i;
 134   unsigned int mode;
 135   double P;
 136   double relerr;
 137   double tmp;
 138
 139   if( t > (n0+n1))
 140     {
 141       GSLEXTRAS_CDF_ERROR("t larger than population size",GSL_EDOM);
 142     }
 143   else if( k >= n0 || k >= t)
 144     {
 145       P = 0.0;
 146     }
 147   else if (k < 0.0)
 148     {
 149       P = 1.0;
 150     }
 151   else
 152     {
 153       P = 0.0;
 154       mode = (int) t*n0 / (n0+n1);
 155       relerr = 1.0;
 156
 157       if(k < mode)
 158         {
 159           i = mode;
 160           while( i <= t && relerr > GSL_DBL_EPSILON && P < 1.0)
 161             {
 162               tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
 163               P += tmp;
 164               relerr = tmp / P;
 165               i++;
 166             }
 167           i = mode - 1;
 168           relerr = 1.0;
 169           while ( i > k && relerr > GSL_DBL_EPSILON && P < 1.0)
 170             {
 171               tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
 172               P += tmp;
 173               relerr = tmp / P;
 174               i--;
 175             }
 176         }
 177       else
 178         {
 179           i = k+1;
 180           while(i <= t && relerr > GSL_DBL_EPSILON && P < 1.0)
 181             {
 182               tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
 183               P += tmp;
 184               relerr = tmp / P;
 185               i++;
 186             }
 187         }
 188       /*
 189        * Hack to get rid of a pesky error when the sum
 190        * gets slightly above 1.0.
 191        */
 192       P = GSL_MIN_DBL(P, 1.0);
 193     }
 194   return P;
 195 }