annotate ezBAMQC/src/htslib/kfunc.c @ 4:50a9d8992e65

Uploaded
author cshl-bsr
date Tue, 29 Mar 2016 15:33:36 -0400
parents dfa3745e5fd8
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
1 /* The MIT License
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
2
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
3 Copyright (C) 2010, 2013 Genome Research Ltd.
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
4 Copyright (C) 2011 Attractive Chaos <attractor@live.co.uk>
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
5
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
6 Permission is hereby granted, free of charge, to any person obtaining
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
7 a copy of this software and associated documentation files (the
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
8 "Software"), to deal in the Software without restriction, including
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
9 without limitation the rights to use, copy, modify, merge, publish,
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
10 distribute, sublicense, and/or sell copies of the Software, and to
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
11 permit persons to whom the Software is furnished to do so, subject to
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
12 the following conditions:
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
13
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
14 The above copyright notice and this permission notice shall be
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
15 included in all copies or substantial portions of the Software.
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
16
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
17 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
18 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
19 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
20 NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
21 BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
22 ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
23 CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
24 SOFTWARE.
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
25 */
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
26
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
27 #include <math.h>
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
28 #include <stdlib.h>
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
29 #include "htslib/kfunc.h"
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
30
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
31 /* Log gamma function
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
32 * \log{\Gamma(z)}
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
33 * AS245, 2nd algorithm, http://lib.stat.cmu.edu/apstat/245
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
34 */
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
35 double kf_lgamma(double z)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
36 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
37 double x = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
38 x += 0.1659470187408462e-06 / (z+7);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
39 x += 0.9934937113930748e-05 / (z+6);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
40 x -= 0.1385710331296526 / (z+5);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
41 x += 12.50734324009056 / (z+4);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
42 x -= 176.6150291498386 / (z+3);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
43 x += 771.3234287757674 / (z+2);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
44 x -= 1259.139216722289 / (z+1);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
45 x += 676.5203681218835 / z;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
46 x += 0.9999999999995183;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
47 return log(x) - 5.58106146679532777 - z + (z-0.5) * log(z+6.5);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
48 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
49
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
50 /* complementary error function
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
51 * \frac{2}{\sqrt{\pi}} \int_x^{\infty} e^{-t^2} dt
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
52 * AS66, 2nd algorithm, http://lib.stat.cmu.edu/apstat/66
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
53 */
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
54 double kf_erfc(double x)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
55 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
56 const double p0 = 220.2068679123761;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
57 const double p1 = 221.2135961699311;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
58 const double p2 = 112.0792914978709;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
59 const double p3 = 33.912866078383;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
60 const double p4 = 6.37396220353165;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
61 const double p5 = .7003830644436881;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
62 const double p6 = .03526249659989109;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
63 const double q0 = 440.4137358247522;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
64 const double q1 = 793.8265125199484;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
65 const double q2 = 637.3336333788311;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
66 const double q3 = 296.5642487796737;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
67 const double q4 = 86.78073220294608;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
68 const double q5 = 16.06417757920695;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
69 const double q6 = 1.755667163182642;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
70 const double q7 = .08838834764831844;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
71 double expntl, z, p;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
72 z = fabs(x) * M_SQRT2;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
73 if (z > 37.) return x > 0.? 0. : 2.;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
74 expntl = exp(z * z * - .5);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
75 if (z < 10. / M_SQRT2) // for small z
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
76 p = expntl * ((((((p6 * z + p5) * z + p4) * z + p3) * z + p2) * z + p1) * z + p0)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
77 / (((((((q7 * z + q6) * z + q5) * z + q4) * z + q3) * z + q2) * z + q1) * z + q0);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
78 else p = expntl / 2.506628274631001 / (z + 1. / (z + 2. / (z + 3. / (z + 4. / (z + .65)))));
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
79 return x > 0.? 2. * p : 2. * (1. - p);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
80 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
81
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
82 /* The following computes regularized incomplete gamma functions.
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
83 * Formulas are taken from Wiki, with additional input from Numerical
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
84 * Recipes in C (for modified Lentz's algorithm) and AS245
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
85 * (http://lib.stat.cmu.edu/apstat/245).
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
86 *
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
87 * A good online calculator is available at:
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
88 *
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
89 * http://www.danielsoper.com/statcalc/calc23.aspx
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
90 *
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
91 * It calculates upper incomplete gamma function, which equals
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
92 * kf_gammaq(s,z)*tgamma(s).
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
93 */
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
94
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
95 #define KF_GAMMA_EPS 1e-14
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
96 #define KF_TINY 1e-290
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
97
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
98 // regularized lower incomplete gamma function, by series expansion
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
99 static double _kf_gammap(double s, double z)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
100 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
101 double sum, x;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
102 int k;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
103 for (k = 1, sum = x = 1.; k < 100; ++k) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
104 sum += (x *= z / (s + k));
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
105 if (x / sum < KF_GAMMA_EPS) break;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
106 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
107 return exp(s * log(z) - z - kf_lgamma(s + 1.) + log(sum));
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
108 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
109 // regularized upper incomplete gamma function, by continued fraction
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
110 static double _kf_gammaq(double s, double z)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
111 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
112 int j;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
113 double C, D, f;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
114 f = 1. + z - s; C = f; D = 0.;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
115 // Modified Lentz's algorithm for computing continued fraction
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
116 // See Numerical Recipes in C, 2nd edition, section 5.2
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
117 for (j = 1; j < 100; ++j) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
118 double a = j * (s - j), b = (j<<1) + 1 + z - s, d;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
119 D = b + a * D;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
120 if (D < KF_TINY) D = KF_TINY;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
121 C = b + a / C;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
122 if (C < KF_TINY) C = KF_TINY;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
123 D = 1. / D;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
124 d = C * D;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
125 f *= d;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
126 if (fabs(d - 1.) < KF_GAMMA_EPS) break;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
127 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
128 return exp(s * log(z) - z - kf_lgamma(s) - log(f));
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
129 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
130
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
131 double kf_gammap(double s, double z)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
132 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
133 return z <= 1. || z < s? _kf_gammap(s, z) : 1. - _kf_gammaq(s, z);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
134 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
135
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
136 double kf_gammaq(double s, double z)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
137 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
138 return z <= 1. || z < s? 1. - _kf_gammap(s, z) : _kf_gammaq(s, z);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
139 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
140
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
141 /* Regularized incomplete beta function. The method is taken from
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
142 * Numerical Recipe in C, 2nd edition, section 6.4. The following web
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
143 * page calculates the incomplete beta function, which equals
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
144 * kf_betai(a,b,x) * gamma(a) * gamma(b) / gamma(a+b):
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
145 *
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
146 * http://www.danielsoper.com/statcalc/calc36.aspx
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
147 */
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
148 static double kf_betai_aux(double a, double b, double x)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
149 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
150 double C, D, f;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
151 int j;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
152 if (x == 0.) return 0.;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
153 if (x == 1.) return 1.;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
154 f = 1.; C = f; D = 0.;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
155 // Modified Lentz's algorithm for computing continued fraction
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
156 for (j = 1; j < 200; ++j) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
157 double aa, d;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
158 int m = j>>1;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
159 aa = (j&1)? -(a + m) * (a + b + m) * x / ((a + 2*m) * (a + 2*m + 1))
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
160 : m * (b - m) * x / ((a + 2*m - 1) * (a + 2*m));
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
161 D = 1. + aa * D;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
162 if (D < KF_TINY) D = KF_TINY;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
163 C = 1. + aa / C;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
164 if (C < KF_TINY) C = KF_TINY;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
165 D = 1. / D;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
166 d = C * D;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
167 f *= d;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
168 if (fabs(d - 1.) < KF_GAMMA_EPS) break;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
169 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
170 return exp(kf_lgamma(a+b) - kf_lgamma(a) - kf_lgamma(b) + a * log(x) + b * log(1.-x)) / a / f;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
171 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
172 double kf_betai(double a, double b, double x)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
173 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
174 return x < (a + 1.) / (a + b + 2.)? kf_betai_aux(a, b, x) : 1. - kf_betai_aux(b, a, 1. - x);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
175 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
176
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
177 #ifdef KF_MAIN
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
178 #include <stdio.h>
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
179 int main(int argc, char *argv[])
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
180 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
181 double x = 5.5, y = 3;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
182 double a, b;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
183 printf("erfc(%lg): %lg, %lg\n", x, erfc(x), kf_erfc(x));
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
184 printf("upper-gamma(%lg,%lg): %lg\n", x, y, kf_gammaq(y, x)*tgamma(y));
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
185 a = 2; b = 2; x = 0.5;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
186 printf("incomplete-beta(%lg,%lg,%lg): %lg\n", a, b, x, kf_betai(a, b, x) / exp(kf_lgamma(a+b) - kf_lgamma(a) - kf_lgamma(b)));
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
187 return 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
188 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
189 #endif
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
190
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
191
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
192 // log\binom{n}{k}
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
193 static double lbinom(int n, int k)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
194 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
195 if (k == 0 || n == k) return 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
196 return lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
197 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
198
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
199 // n11 n12 | n1_
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
200 // n21 n22 | n2_
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
201 //-----------+----
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
202 // n_1 n_2 | n
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
203
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
204 // hypergeometric distribution
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
205 static double hypergeo(int n11, int n1_, int n_1, int n)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
206 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
207 return exp(lbinom(n1_, n11) + lbinom(n-n1_, n_1-n11) - lbinom(n, n_1));
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
208 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
209
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
210 typedef struct {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
211 int n11, n1_, n_1, n;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
212 double p;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
213 } hgacc_t;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
214
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
215 // incremental version of hypergenometric distribution
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
216 static double hypergeo_acc(int n11, int n1_, int n_1, int n, hgacc_t *aux)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
217 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
218 if (n1_ || n_1 || n) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
219 aux->n11 = n11; aux->n1_ = n1_; aux->n_1 = n_1; aux->n = n;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
220 } else { // then only n11 changed; the rest fixed
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
221 if (n11%11 && n11 + aux->n - aux->n1_ - aux->n_1) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
222 if (n11 == aux->n11 + 1) { // incremental
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
223 aux->p *= (double)(aux->n1_ - aux->n11) / n11
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
224 * (aux->n_1 - aux->n11) / (n11 + aux->n - aux->n1_ - aux->n_1);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
225 aux->n11 = n11;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
226 return aux->p;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
227 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
228 if (n11 == aux->n11 - 1) { // incremental
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
229 aux->p *= (double)aux->n11 / (aux->n1_ - n11)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
230 * (aux->n11 + aux->n - aux->n1_ - aux->n_1) / (aux->n_1 - n11);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
231 aux->n11 = n11;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
232 return aux->p;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
233 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
234 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
235 aux->n11 = n11;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
236 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
237 aux->p = hypergeo(aux->n11, aux->n1_, aux->n_1, aux->n);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
238 return aux->p;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
239 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
240
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
241 double kt_fisher_exact(int n11, int n12, int n21, int n22, double *_left, double *_right, double *two)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
242 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
243 int i, j, max, min;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
244 double p, q, left, right;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
245 hgacc_t aux;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
246 int n1_, n_1, n;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
247
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
248 n1_ = n11 + n12; n_1 = n11 + n21; n = n11 + n12 + n21 + n22; // calculate n1_, n_1 and n
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
249 max = (n_1 < n1_) ? n_1 : n1_; // max n11, for right tail
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
250 min = n1_ + n_1 - n; // not sure why n11-n22 is used instead of min(n_1,n1_)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
251 if (min < 0) min = 0; // min n11, for left tail
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
252 *two = *_left = *_right = 1.;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
253 if (min == max) return 1.; // no need to do test
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
254 q = hypergeo_acc(n11, n1_, n_1, n, &aux); // the probability of the current table
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
255 // left tail
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
256 p = hypergeo_acc(min, 0, 0, 0, &aux);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
257 for (left = 0., i = min + 1; p < 0.99999999 * q && i<=max; ++i) // loop until underflow
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
258 left += p, p = hypergeo_acc(i, 0, 0, 0, &aux);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
259 --i;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
260 if (p < 1.00000001 * q) left += p;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
261 else --i;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
262 // right tail
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
263 p = hypergeo_acc(max, 0, 0, 0, &aux);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
264 for (right = 0., j = max - 1; p < 0.99999999 * q && j>=0; --j) // loop until underflow
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
265 right += p, p = hypergeo_acc(j, 0, 0, 0, &aux);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
266 ++j;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
267 if (p < 1.00000001 * q) right += p;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
268 else ++j;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
269 // two-tail
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
270 *two = left + right;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
271 if (*two > 1.) *two = 1.;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
272 // adjust left and right
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
273 if (abs(i - n11) < abs(j - n11)) right = 1. - left + q;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
274 else left = 1.0 - right + q;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
275 *_left = left; *_right = right;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
276 return q;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
277 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
278
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
279
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
280