1 /*
2 * AUTHOR
3 * Catherine Loader, catherine@research.bell-labs.com.
4 * October 23, 2000.
5 *
6 * Merge in to R:
7 * Copyright (C) 2000, 2001 The R Core Team
8 *
9 * This program is free software; you can redistribute it and/or modify
10 * it under the terms of the GNU General Public License as published by
11 * the Free Software Foundation; either version 2 of the License, or
12 * (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, a copy is available at
21 * http://www.r-project.org/Licenses/
22 *
23 *
24 * DESCRIPTION
25 *
26 * Given a sequence of r successes and b failures, we sample n (\le b+r)
27 * items without replacement. The hypergeometric probability is the
28 * probability of x successes:
29 *
30 * choose(r, x) * choose(b, n-x)
31 * p(x; r,b,n) = ----------------------------- =
32 * choose(r+b, n)
33 *
34 * dbinom(x,r,p) * dbinom(n-x,b,p)
35 * = --------------------------------
36 * dbinom(n,r+b,p)
37 *
38 * for any p. For numerical stability, we take p=n/(r+b); with this choice,
39 * the denominator is not exponentially small.
40 */
41
42 #include "nmath.h"
43 #include "dpq.h"
44
dhyper(double x,double r,double b,double n,int give_log)45 double dhyper(double x, double r, double b, double n, int give_log)
46 {
47 double p, q, p1, p2, p3;
48
49 #ifdef IEEE_754
50 if (ISNAN(x) || ISNAN(r) || ISNAN(b) || ISNAN(n))
51 return x + r + b + n;
52 #endif
53
54 if (R_D_negInonint(r) || R_D_negInonint(b) || R_D_negInonint(n) || n > r+b)
55 ML_ERR_return_NAN;
56 if (R_D_negInonint(x))
57 return(R_D__0);
58
59 x = R_D_forceint(x);
60 r = R_D_forceint(r);
61 b = R_D_forceint(b);
62 n = R_D_forceint(n);
63
64 if (n < x || r < x || n - x > b) return(R_D__0);
65 if (n == 0) return((x == 0) ? R_D__1 : R_D__0);
66
67 p = ((double)n)/((double)(r+b));
68 q = ((double)(r+b-n))/((double)(r+b));
69
70 p1 = dbinom_raw(x, r, p,q,give_log);
71 p2 = dbinom_raw(n-x,b, p,q,give_log);
72 p3 = dbinom_raw(n,r+b, p,q,give_log);
73
74 return( (give_log) ? p1 + p2 - p3 : p1*p2/p3 );
75 }
76