1 /* randist/hyperg.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include "gsl__config.h"
21 #include <math.h>
22 #include "gsl_rng.h"
23 #include "gsl_randist.h"
24 #include "gsl_sf_gamma.h"
25
26 /* The hypergeometric distribution has the form,
27
28 prob(k) = choose(n1,t) choose(n2, t-k) / choose(n1+n2,t)
29
30 where choose(a,b) = a!/(b!(a-b)!)
31
32 n1 + n2 is the total population (tagged plus untagged)
33 n1 is the tagged population
34 t is the number of samples taken (without replacement)
35 k is the number of tagged samples found
36
37 */
38
39 unsigned int
gsl_ran_hypergeometric(const gsl_rng * r,unsigned int n1,unsigned int n2,unsigned int t)40 gsl_ran_hypergeometric (const gsl_rng * r, unsigned int n1, unsigned int n2,
41 unsigned int t)
42 {
43 const unsigned int n = n1 + n2;
44
45 unsigned int i = 0;
46 unsigned int a = n1;
47 unsigned int b = n1 + n2;
48 unsigned int k = 0;
49
50 if (t > n)
51 {
52 t = n ;
53 }
54
55 if (t < n / 2)
56 {
57 for (i = 0 ; i < t ; i++)
58 {
59 double u = gsl_rng_uniform(r) ;
60
61 if (b * u < a)
62 {
63 k++ ;
64 if (k == n1)
65 return k ;
66 a-- ;
67 }
68 b-- ;
69 }
70 return k;
71 }
72 else
73 {
74 for (i = 0 ; i < n - t ; i++)
75 {
76 double u = gsl_rng_uniform(r) ;
77
78 if (b * u < a)
79 {
80 k++ ;
81 if (k == n1)
82 return n1 - k ;
83 a-- ;
84 }
85 b-- ;
86 }
87 return n1 - k;
88 }
89
90
91 }
92
93 double
gsl_ran_hypergeometric_pdf(const unsigned int k,const unsigned int n1,const unsigned int n2,unsigned int t)94 gsl_ran_hypergeometric_pdf (const unsigned int k,
95 const unsigned int n1,
96 const unsigned int n2,
97 unsigned int t)
98 {
99 if (t > n1 + n2)
100 {
101 t = n1 + n2 ;
102 }
103
104 if (k > n1 || k > t)
105 {
106 return 0 ;
107 }
108 else if (t > n2 && k + n2 < t )
109 {
110 return 0 ;
111 }
112 else
113 {
114 double p;
115
116 double c1 = gsl_sf_lnchoose(n1,k);
117 double c2 = gsl_sf_lnchoose(n2,t-k);
118 double c3 = gsl_sf_lnchoose(n1+n2,t);
119
120 p = exp(c1 + c2 - c3) ;
121
122 return p;
123 }
124 }
125