|
0
|
1 #include <math.h>
|
|
|
2 #include <stdlib.h>
|
|
|
3
|
|
|
4 /* This program is implemented with ideas from this web page:
|
|
|
5 *
|
|
|
6 * http://www.langsrud.com/fisher.htm
|
|
|
7 */
|
|
|
8
|
|
|
9 // log\binom{n}{k}
|
|
|
10 static double lbinom(int n, int k)
|
|
|
11 {
|
|
|
12 if (k == 0 || n == k) return 0;
|
|
|
13 return lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1);
|
|
|
14 }
|
|
|
15
|
|
|
16 // n11 n12 | n1_
|
|
|
17 // n21 n22 | n2_
|
|
|
18 //-----------+----
|
|
|
19 // n_1 n_2 | n
|
|
|
20
|
|
|
21 // hypergeometric distribution
|
|
|
22 static double hypergeo(int n11, int n1_, int n_1, int n)
|
|
|
23 {
|
|
|
24 return exp(lbinom(n1_, n11) + lbinom(n-n1_, n_1-n11) - lbinom(n, n_1));
|
|
|
25 }
|
|
|
26
|
|
|
27 typedef struct {
|
|
|
28 int n11, n1_, n_1, n;
|
|
|
29 double p;
|
|
|
30 } hgacc_t;
|
|
|
31
|
|
|
32 // incremental version of hypergenometric distribution
|
|
|
33 static double hypergeo_acc(int n11, int n1_, int n_1, int n, hgacc_t *aux)
|
|
|
34 {
|
|
|
35 if (n1_ || n_1 || n) {
|
|
|
36 aux->n11 = n11; aux->n1_ = n1_; aux->n_1 = n_1; aux->n = n;
|
|
|
37 } else { // then only n11 changed; the rest fixed
|
|
|
38 if (n11%11 && n11 + aux->n - aux->n1_ - aux->n_1) {
|
|
|
39 if (n11 == aux->n11 + 1) { // incremental
|
|
|
40 aux->p *= (double)(aux->n1_ - aux->n11) / n11
|
|
|
41 * (aux->n_1 - aux->n11) / (n11 + aux->n - aux->n1_ - aux->n_1);
|
|
|
42 aux->n11 = n11;
|
|
|
43 return aux->p;
|
|
|
44 }
|
|
|
45 if (n11 == aux->n11 - 1) { // incremental
|
|
|
46 aux->p *= (double)aux->n11 / (aux->n1_ - n11)
|
|
|
47 * (aux->n11 + aux->n - aux->n1_ - aux->n_1) / (aux->n_1 - n11);
|
|
|
48 aux->n11 = n11;
|
|
|
49 return aux->p;
|
|
|
50 }
|
|
|
51 }
|
|
|
52 aux->n11 = n11;
|
|
|
53 }
|
|
|
54 aux->p = hypergeo(aux->n11, aux->n1_, aux->n_1, aux->n);
|
|
|
55 return aux->p;
|
|
|
56 }
|
|
|
57
|
|
|
58 double kt_fisher_exact(int n11, int n12, int n21, int n22, double *_left, double *_right, double *two)
|
|
|
59 {
|
|
|
60 int i, j, max, min;
|
|
|
61 double p, q, left, right;
|
|
|
62 hgacc_t aux;
|
|
|
63 int n1_, n_1, n;
|
|
|
64
|
|
|
65 n1_ = n11 + n12; n_1 = n11 + n21; n = n11 + n12 + n21 + n22; // calculate n1_, n_1 and n
|
|
|
66 max = (n_1 < n1_) ? n_1 : n1_; // max n11, for right tail
|
|
|
67 min = n1_ + n_1 - n;
|
|
|
68 if (min < 0) min = 0; // min n11, for left tail
|
|
|
69 *two = *_left = *_right = 1.;
|
|
|
70 if (min == max) return 1.; // no need to do test
|
|
|
71 q = hypergeo_acc(n11, n1_, n_1, n, &aux); // the probability of the current table
|
|
|
72 // left tail
|
|
|
73 p = hypergeo_acc(min, 0, 0, 0, &aux);
|
|
|
74 for (left = 0., i = min + 1; p < 0.99999999 * q; ++i) // loop until underflow
|
|
|
75 left += p, p = hypergeo_acc(i, 0, 0, 0, &aux);
|
|
|
76 --i;
|
|
|
77 if (p < 1.00000001 * q) left += p;
|
|
|
78 else --i;
|
|
|
79 // right tail
|
|
|
80 p = hypergeo_acc(max, 0, 0, 0, &aux);
|
|
|
81 for (right = 0., j = max - 1; p < 0.99999999 * q; --j) // loop until underflow
|
|
|
82 right += p, p = hypergeo_acc(j, 0, 0, 0, &aux);
|
|
|
83 ++j;
|
|
|
84 if (p < 1.00000001 * q) right += p;
|
|
|
85 else ++j;
|
|
|
86 // two-tail
|
|
|
87 *two = left + right;
|
|
|
88 if (*two > 1.) *two = 1.;
|
|
|
89 // adjust left and right
|
|
|
90 if (abs(i - n11) < abs(j - n11)) right = 1. - left + q;
|
|
|
91 else left = 1.0 - right + q;
|
|
|
92 *_left = left; *_right = right;
|
|
|
93 return q;
|
|
|
94 }
|
|
|
95
|
|
|
96 #ifdef FET_MAIN
|
|
|
97 #include <stdio.h>
|
|
|
98
|
|
|
99 int main(int argc, char *argv[])
|
|
|
100 {
|
|
|
101 char id[1024];
|
|
|
102 int n11, n12, n21, n22;
|
|
|
103 double left, right, twotail, prob;
|
|
|
104
|
|
|
105 while (scanf("%s%d%d%d%d", id, &n11, &n12, &n21, &n22) == 5) {
|
|
|
106 prob = kt_fisher_exact(n11, n12, n21, n22, &left, &right, &twotail);
|
|
|
107 printf("%s\t%d\t%d\t%d\t%d\t%.6g\t%.6g\t%.6g\t%.6g\n", id, n11, n12, n21, n22,
|
|
|
108 prob, left, right, twotail);
|
|
|
109 }
|
|
|
110 return 0;
|
|
|
111 }
|
|
|
112 #endif
|
