source: orange-bioinformatics/orangecontrib/bio/obiProb.py @ 1890:043f2b1864cd

Revision 1890:043f2b1864cd, 5.6 KB checked in by markotoplak, 6 months ago (diff)

Numberically stable p-values for obiProb.Binomial.

Line 
1import math   
2
3def _lngamma(z):
4    x = 0
5    x += 0.1659470187408462e-06 / (z + 7)
6    x += 0.9934937113930748e-05 / (z + 6)
7    x -= 0.1385710331296526 / (z + 5)
8    x += 12.50734324009056 / (z + 4)
9    x -= 176.6150291498386 / (z + 3)
10    x += 771.3234287757674 / (z + 2)
11    x -= 1259.139216722289 / (z + 1)
12    x += 676.5203681218835 / (z)
13    x += 0.9999999999995183
14   
15    return math.log(x) - 5.58106146679532777 - z + (z - 0.5) * math.log(z + 6.5)
16       
17class LogBin(object):
18    _max = 2
19    _lookup = [0.0, 0.0]
20    _max_factorial = 1
21    def __init__(self, max=1000):
22        self._extend(max)
23
24    @staticmethod
25    def _extend(max):
26        if max <= LogBin._max:
27            return
28        for i in range(LogBin._max, max):
29            if i > 1000: ## an arbitrary cuttof
30                LogBin._lookup.append(LogBin._logfactorial(i))
31            else:
32                LogBin._max_factorial *= i
33                LogBin._lookup.append(math.log(LogBin._max_factorial))
34        LogBin._max = max
35
36    def _logbin(self, n, k):
37        if n >= self._max:
38            self._extend(n + 100)
39        if k < n and k >= 0:
40            return self._lookup[n] - self._lookup[n - k] - self._lookup[k]
41        else:
42            return 0.0
43
44    @staticmethod
45    def _logfactorial(n):
46        if (n <= 1):
47            return 0.0
48        else:
49            return _lngamma(n + 1)
50
51class Binomial(LogBin):
52
53    def __call__(self, k, N, m, n):
54        p = 1.0 * m / N
55        if p == 0.0:
56            if k == 0:
57                return 1.0
58            else:
59                return 0.0
60        elif p == 1.0:
61            if n == k:
62                return 1.0
63            else:
64                return 0.0
65        try:
66            return min(math.exp(self._logbin(n, k) + k * math.log(p) + (n - k) * math.log(1.0 - p)), 1.0)
67        except (OverflowError, ValueError), er:
68            print k, N, m, n
69            raise
70##        return math.exp(self._logbin(n, k) + math.log((p**k) * (1.0 - p)**(n - k)))
71
72    def p_value(self, k, N, m, n):
73        if n - k + 1 <= k:
74            #starting from k gives the shorter list of values
75            return sum(self.__call__(i, N, m, n) for i in range(k, n+1))
76        else:
77            value = 1.0 - sum(self.__call__(i, N, m, n) for i in range(k))
78            #if the value is small it is probably inexact due to the limited
79            #precision of floats, as for example  (1-(1-1e-20)) -> 0
80            #if so, compute the result without substraction
81            if value < 1e-3: #arbitary threshold
82                #print "INEXACT", value, sum(self.__call__(i, N, m, n) for i in range(k, n+1))
83                return sum(self.__call__(i, N, m, n) for i in range(k, n+1))
84            else:
85                return value
86
87class Hypergeometric(LogBin):
88
89    def __call__(self, k, N, m, n):
90        if k < max(0, n + m - N) or k > min(n, m):
91            return 0.0
92        try:
93            return min(math.exp(self._logbin(m, k) + self._logbin(N - m, n - k) - self._logbin(N, n)), 1.0)
94        except (OverflowError, ValueError), er:
95            print k, N, m, n
96            raise
97
98    def p_value(self, k, N, m, n):
99        #From wolfram alpha:
100        #1-CDF[HypergeometricDistribution[80, 50, 2000], 20-1] = 1.75695e-16
101        #1-CDF[HypergeometricDistribution[80, 50, 2000], 40-1] = 9.92008e-52
102
103        if min(n,m) - k + 1 <= k:
104            #starting from k gives the shorter list of values
105            return sum(self.__call__(i, N, m, n) for i in range(k, min(n,m)+1))
106        else:
107            value = 1.0 - sum(self.__call__(i, N, m, n) for i in (range(k)))
108            #if the value is small it is probably inexact due to the limited
109            #precision of floats, as for example  (1-(1-1e-20)) -> 0
110            #if so, compute the result without substraction
111            if value < 1e-3: #arbitary threshold
112                #print "INEXACT", value, sum(self.__call__(i, N, m, n) for i in range(k, min(n,m)+1))
113                return sum(self.__call__(i, N, m, n) for i in range(k, min(n,m)+1))
114            else:
115                return value
116
117## to speed-up FDR, calculate ahead sum([1/i for i in range(1, m+1)]), for m in [1,100000]. For higher values of m use an approximation, with error less or equal to 4.99999157277e-006. (sum([1/i for i in range(1, m+1)])  ~ log(m) + 0.5772..., 0.5572 is an Euler-Mascheroni constant)
118c = [1.0]
119for m in range(2, 100000):
120    c.append( c[-1] + 1.0/m)
121
122def is_sorted(l):
123    return all(l[i] <= l[i+1] for i in xrange(len(l)-1))
124
125def FDR(p_values, dependent=False, m=None, ordered=False):
126    """
127    If the user is sure that pvalues as already sorted nondescendingly
128    setting ordered=True will make the computation faster.
129    """
130
131    if not ordered:
132        ordered = is_sorted(p_values)
133
134    if not ordered:
135        joined = [ (v,i) for i,v in enumerate(p_values) ]
136        joined.sort()
137        p_values = [ p[0] for p in joined ]
138        indices = [ p[1] for p in joined ]
139
140    if not m:
141        m = len(p_values)
142    if m <= 0 or not p_values:
143        return []
144
145    if dependent: # correct q for dependent tests
146        k = c[m-1] if m <= len(c) else math.log(m) + 0.57721566490153286060651209008240243104215933593992
147        m = m * k
148
149    tmp_fdrs = [p*m/(i+1.0) for (i, p) in enumerate(p_values)]
150    fdrs = []
151    cmin = tmp_fdrs[-1]
152    for f in reversed(tmp_fdrs):
153        cmin = min(f, cmin)
154        fdrs.append( cmin)
155    fdrs.reverse()
156
157    if not ordered:
158        new = [ None ] * len(fdrs)
159        for v,i in zip(fdrs, indices):
160            new[i] = v
161        fdrs = new
162
163    return fdrs
164
165def Bonferroni(p_values, m=None):
166    if not m:
167        m = len(p_values)
168    if m == 0:
169        return []
170    m = float(m)
171    return [p/m for p in p_values]
Note: See TracBrowser for help on using the repository browser.