Fix PostScript output of multiple charts.
[pspp] / 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 <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 }