| 42303 |
ripley |
1 |
/*
|
|
|
2 |
* AUTHOR
|
|
|
3 |
* Catherine Loader, catherine@research.bell-labs.com.
|
|
|
4 |
* October 23, 2000.
|
|
|
5 |
*
|
|
|
6 |
* Merge in to R:
|
| 64668 |
ripley |
7 |
* Copyright (C) 2000-2014 The R Core Team
|
| 42303 |
ripley |
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
|
| 68947 |
ripley |
21 |
* https://www.R-project.org/Licenses/
|
| 42303 |
ripley |
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 |
|
|
|
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;
|
| 65344 |
maechler |
56 |
if(x < 0) return(R_D__0);
|
|
|
57 |
R_D_nonint_check(x);// incl warning
|
| 42303 |
ripley |
58 |
|
| 64668 |
ripley |
59 |
x = R_forceint(x);
|
|
|
60 |
r = R_forceint(r);
|
|
|
61 |
b = R_forceint(b);
|
|
|
62 |
n = R_forceint(n);
|
| 42303 |
ripley |
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 |
}
|
The R Project SVN
R