1 /* cdf/hypergeometric.c
3 * Copyright (C) 2004 Free Software Foundation, Inc.
4 * Written by Jason H. Stover.
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.
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.
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.
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.
27 * This algorithm computes Pr( X <= k ) by summing the terms from
28 * the mass function, Pr( X = k ).
32 * T. Wu. An accurate computation of the hypergeometric distribution
33 * function. ACM Transactions on Mathematical Software. Volume 19, number 1,
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.
40 * W. Feller. An Introduction to Probability Theory and Its Applications,
41 * third edition. 1968. Chapter 2, section 6.
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"
55 gslextras_cdf_hypergeometric_P (const unsigned int k,
56 const unsigned int n0,
57 const unsigned int n1,
68 GSLEXTRAS_CDF_ERROR("t larger than population size",GSL_EDOM);
70 else if( k >= n0 || k >= t)
81 mode = (int) t*n0 / (n0+n1);
87 while(i != UINT_MAX && relerr > GSL_DBL_EPSILON && P < 1.0)
89 tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
99 while(i <= k && relerr > GSL_DBL_EPSILON && P < 1.0)
101 tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
108 while( i != UINT_MAX && relerr > GSL_DBL_EPSILON && P < 1.0)
110 tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
117 * Hack to get rid of a pesky error when the sum
118 * gets slightly above 1.0.
120 P = GSL_MIN_DBL (P, 1.0);
129 gslextras_cdf_hypergeometric_Q (const unsigned int k,
130 const unsigned int n0,
131 const unsigned int n1,
132 const unsigned int t)
142 GSLEXTRAS_CDF_ERROR("t larger than population size",GSL_EDOM);
144 else if( k >= n0 || k >= t)
155 mode = (int) t*n0 / (n0+n1);
161 while( i <= t && relerr > GSL_DBL_EPSILON && P < 1.0)
163 tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
170 while ( i > k && relerr > GSL_DBL_EPSILON && P < 1.0)
172 tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
181 while(i <= t && relerr > GSL_DBL_EPSILON && P < 1.0)
183 tmp = gsl_ran_hypergeometric_pdf(i, n0, n1, t);
190 * Hack to get rid of a pesky error when the sum
191 * gets slightly above 1.0.
193 P = GSL_MIN_DBL(P, 1.0);