lib/gsl-extras/hypergeometric.c

   1 /* cdf/hypergeometric.c
   2  *
   3  * Copyright (C) 2004 Free Software Foundation, Inc.
   4  * Written by Jason H. Stover.
   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 2 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  19  */
  20
  21 /*
  22  * Computes the cumulative distribution function for a hypergeometric
  23  * random variable. A hypergeometric random variable X is the number
  24  * of elements of type 0 in a sample of size t, drawn from a population
  25  * of size n1 + n0, in which n1 are of type 1 and n0 are of type 0.
  26  *
  27  * This algorithm computes Pr( X <= k ) by summing the terms from
  28  * the mass function, Pr( X = k ).
  29  *
  30  * References:
  31  *
  32  * T. Wu. An accurate computation of the hypergeometric distribution
  33  * function. ACM Transactions on Mathematical Software. Volume 19, number 1,
  34  * March 1993.
  35  *  This algorithm is not used, since it requires factoring the
  36  *  numerator and denominator, then cancelling. It is more accurate
  37  *  than the algorithm used here, but the cancellation requires more
  38  *  time than the algorithm used here.
  39  *
  40  * W. Feller. An Introduction to Probability Theory and Its Applications,
  41  * third edition. 1968. Chapter 2, section 6.
  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 #include "gsl-extras.h"
  50
  51 /*
  52  * Pr (X <= k)
  53  */
  54 double
  55 gslextras_cdf_hypergeometric_P (const unsigned int k,
  56                                 const unsigned int n0,
  57                                 const unsigned int n1,
  58                                 const unsigned int t)
  59 {
  60   unsigned int i;
  61   unsigned int mode;
  62   double P;
  63   double tmp;
  64   double relerr;
  65
  66   if( t > (n0+n1))
  67     {
  68       GSLEXTRAS_CDF_ERROR("t larger than population size",GSL_EDOM);
  69     }
  70   else if( k >= n0 || k >= t)
  71     {
  72       P = 1.0;
  73     }
  74   else if (k < 0.0)
  75     {
  76       P = 0.0;
  77     }
  78   else
  79     {
  80       P = 0.0;
  81       mode = (int) t*n0 / (n0+n1);
  82       relerr = 1.0;
  83       if( k < mode )
  84         {
  85           i = k;
  86           relerr = 1.0;
  87           while(i != UINT_MAX && relerr > GSL_DBL_EPSILON && P < 1.0)
  88             {
  89               tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
  90               P += tmp;
  91               relerr = tmp / P;
  92               i--;
  93             }
  94         }
  95       else
  96         {
  97           i = mode;
  98           relerr = 1.0;
  99           while(i <= k && relerr > GSL_DBL_EPSILON && P < 1.0)
 100             {
 101               tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
 102               P += tmp;
 103               relerr = tmp / P;
 104               i++;
 105             }
 106           i = mode - 1;
 107           relerr = 1.0;
 108           while( i != UINT_MAX && relerr > GSL_DBL_EPSILON && P < 1.0)
 109             {
 110               tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
 111               P += tmp;
 112               relerr = tmp / P;
 113               i--;
 114             }
 115         }
 116       /*
 117        * Hack to get rid of a pesky error when the sum
 118        * gets slightly above 1.0.
 119        */
 120       P = GSL_MIN_DBL (P, 1.0);
 121     }
 122   return P;
 123 }
 124
 125 /*
 126  * Pr (X > k)
 127  */
 128 double
 129 gslextras_cdf_hypergeometric_Q (const unsigned int k,
 130                                 const unsigned int n0,
 131                                 const unsigned int n1,
 132                                 const unsigned int t)
 133 {
 134   unsigned int i;
 135   unsigned int mode;
 136   double P;
 137   double relerr;
 138   double tmp;
 139
 140   if( t > (n0+n1))
 141     {
 142       GSLEXTRAS_CDF_ERROR("t larger than population size",GSL_EDOM);
 143     }
 144   else if( k >= n0 || k >= t)
 145     {
 146       P = 0.0;
 147     }
 148   else if (k < 0.0)
 149     {
 150       P = 1.0;
 151     }
 152   else
 153     {
 154       P = 0.0;
 155       mode = (int) t*n0 / (n0+n1);
 156       relerr = 1.0;
 157
 158       if(k < mode)
 159         {
 160           i = mode;
 161           while( i <= t && relerr > GSL_DBL_EPSILON && P < 1.0)
 162             {
 163               tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
 164               P += tmp;
 165               relerr = tmp / P;
 166               i++;
 167             }
 168           i = mode - 1;
 169           relerr = 1.0;
 170           while ( i > k && relerr > GSL_DBL_EPSILON && P < 1.0)
 171             {
 172               tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
 173               P += tmp;
 174               relerr = tmp / P;
 175               i--;
 176             }
 177         }
 178       else
 179         {
 180           i = k+1;
 181           while(i <= t && relerr > GSL_DBL_EPSILON && P < 1.0)
 182             {
 183               tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
 184               P += tmp;
 185               relerr = tmp / P;
 186               i++;
 187             }
 188         }
 189       /*
 190        * Hack to get rid of a pesky error when the sum
 191        * gets slightly above 1.0.
 192        */
 193       P = GSL_MIN_DBL(P, 1.0);
 194     }
 195   return P;
 196 }