source: orange/orange/orngContingency.py @ 6538:a5f65d7f0b2c

Revision 6538:a5f65d7f0b2c, 10.5 KB checked in by Mitar <Mitar@…>, 4 years ago (diff)

Made XPM version of the icon 32x32.

Line 
1#
2# Module Orange Contingency Tables
3# --------------------------------
4#
5# CVS Status: $Id$
6#
7# Author: Aleks Jakulin (jakulin@acm.org)
8# (Copyright (C)2004 Aleks Jakulin)
9#
10# Purpose: Management of contingency tables and various heuristics.
11#          For performance reasons, handling of 2 and 3-way tables
12#          is handled separately.
13#
14# Project initiated on 2003/01/30
15#
16# ChangeLog:
17#
18
19
20import numpy, orange, statc
21from random import *
22
23_log2 = 1.0/numpy.log(2)
24_log2e = numpy.log(2)
25
26def Flatten(m):
27    if len(numpy.shape(m)) > 1:
28        v = numpy.ravel(m)
29    else:
30        v = m
31    return v
32
33def Probabilities(m):
34    t = numpy.sum(Flatten(m))
35    if t == 0:
36        return 0
37    else:
38        return m/t
39
40def Entropy(m):
41    v = Flatten(m)
42    pv = Probabilities(v)
43    lv = numpy.log(numpy.clip(pv,1e-6,1.0))*_log2
44    return -numpy.dot(pv,lv)
45
46
47class ContingencyTable3:
48    def InteractionInformation(self):
49        abc = Entropy(self.m)
50        return (-Entropy(self.a)-Entropy(self.b)-Entropy(self.c)+Entropy(self.ab)+Entropy(self.bc)+Entropy(self.ac)-abc)
51
52    def CMI(self):
53        # the three conditional mutual informations
54        abc = Entropy(self.m)
55        c = (-Entropy(self.c)+Entropy(self.bc)+Entropy(self.ac)-abc)
56        b = (-Entropy(self.b)+Entropy(self.ab)+Entropy(self.bc)-abc)
57        a = (-Entropy(self.a)+Entropy(self.ab)+Entropy(self.ac)-abc)
58        return (a,b,c)
59
60    def JaccardInteraction(self):
61        abc = Entropy(self.m)
62        return (-Entropy(self.a)-Entropy(self.b)-Entropy(self.c)+Entropy(self.ab)+Entropy(self.bc)+Entropy(self.ac)-abc)/abc
63
64    def NormDivergence(self):
65        try:
66            return self.kirknorm
67        except:
68            s = 0.0
69            for x in range(len(self.values[0])):
70                for y in range(len(self.values[1])):
71                    for z in range(len(self.values[2])):
72                        s += self.Divergence(x,y,z)[1]
73            self.kirknorm = s
74            return s
75
76    def IPF(self,tolerance=1e-6,maxiterations=100):
77        d = numpy.shape(self.m)
78        p = self.pm
79        px = self.pab
80        py = self.pac
81        pz = self.pbc
82        pxx = self.pa
83        pyy = self.pb
84        pzz = self.pc
85        v = 1.0/(d[0]*d[1]*d[2])
86        model = numpy.ones(d,numpy.float)*v
87
88        iterations = 0
89        diff = 1e30
90        pdiv = 1e50
91        while iterations < maxiterations and abs(diff) > tolerance:
92            # FIT
93            for c in range(3):
94                if c == 0:
95                    mx = Probabilities(numpy.sum(model,axis=2))
96                elif c == 1:
97                    my = Probabilities(numpy.sum(model,axis=1))
98                else:
99                    mz = Probabilities(numpy.sum(model,axis=0))
100                for x in xrange(d[0]):
101                    for y in xrange(d[1]):
102                        for z in xrange(d[2]):
103                            if c == 0:
104                                model[x,y,z] *= px[x,y]/max(mx[x,y],1e-16)
105                            elif c == 1:
106                                model[x,y,z] *= py[x,z]/max(my[x,z],1e-16)
107                            else:
108                                model[x,y,z] *= pz[y,z]/max(mz[y,z],1e-16)
109            # EVALUATE
110            div = 0.0
111            for x in xrange(d[0]):
112                for y in xrange(d[1]):
113                    for z in xrange(d[2]):
114                        if p[x,y,z]>0:
115                            div += p[x,y,z]*numpy.log(p[x,y,z]/model[x,y,z])
116            #print iterations, Entropy(model), div
117            iterations += 1
118            diff = pdiv-div
119            pdiv = div
120        self.gis = model
121        return div*_log2
122
123    def KSA(self):
124        d = numpy.shape(self.m)
125        self.kirk = numpy.ones(d,numpy.float)
126
127        # normalize kirkwood approximation
128        sumx = 0.0
129        for x in xrange(d[0]):
130            for y in xrange(d[1]):
131                for z in xrange(d[2]):
132                    pkirkwood  = self.pab[x,y]*self.pbc[y,z]*self.pac[x,z]
133                    q = self.pa[x]*self.pb[y]*self.pc[z]
134                    if q > 0:
135                        pkirkwood /= self.pa[x]*self.pb[y]*self.pc[z]
136                    else:
137                        pkirkwood = 0.0
138                    sumx += pkirkwood
139                    self.kirk[x,y,z] = pkirkwood
140        self.kirk /= sumx
141
142        div = 0.0
143        for x in xrange(d[0]):
144            for y in xrange(d[1]):
145                for z in xrange(d[2]):
146                    if self.pm[x,y,z]>0:
147                        div += self.pm[x,y,z]*numpy.log(self.pm[x,y,z]/self.kirk[x,y,z])
148        return (div*_log2,sumx)
149
150    def Divergence(self,x,y,z):
151        ptrue = self.pm[x,y,z]
152        pkirkwood  = self.pab[x,y]*self.pbc[y,z]*self.pac[x,z]
153        pkirkwood /= self.pa[x]*self.pb[y]*self.pc[z]
154        if ptrue > 1e-6:
155            div = numpy.log(ptrue/pkirkwood)
156        else:
157            div = 0.0
158        return (ptrue,pkirkwood,_log2*div)
159
160    def NDivergence(self,x,y,z):
161        norm = 1.0/self.NormDivergence()
162        ptrue = self.pm[x,y,z]
163        pkirkwood  = self.pab[x,y]*self.pbc[y,z]*self.pac[x,z]
164        pkirkwood /= self.pa[x]*self.pb[y]*self.pc[z]
165        if ptrue > 1e-6:
166            div = numpy.log(ptrue/(pkirkwood*norm))
167        else:
168            div = 0.0
169        return (ptrue,pkirkwood*norm,_log2*div)
170
171    def __init__(self, m, names, values):
172        self.names  = names
173        self.values = values
174        m = numpy.array(m,numpy.float)
175        self.m = m
176        self.bc = numpy.sum(m,axis=0)
177        self.ac = numpy.sum(m,axis=1)
178        self.ab = numpy.sum(m,axis=2)
179
180        self.a = numpy.sum(self.ab,axis=1)
181        self.b = numpy.sum(self.ab,axis=0)
182        self.c = numpy.sum(self.ac,axis=0)
183        self.total = numpy.sum(self.a)
184
185        self.pm = Probabilities(self.m)
186        self.pab = Probabilities(self.ab)
187        self.pbc = Probabilities(self.bc)
188        self.pac = Probabilities(self.ac)
189        self.pa = Probabilities(self.a)
190        self.pb = Probabilities(self.b)
191        self.pc = Probabilities(self.c)
192        dof = 0
193        (ni,nj,nk) = numpy.shape(self.m)
194        for ii in xrange(ni):
195            for jj in xrange(nj):
196                for kk in xrange(nk):
197                    if self.m[ii,jj,kk] > 0:
198                        dof += 1
199        self.dof = dof-1
200        return
201
202class ContingencyTable2:
203    def InteractionInformation(self):
204        return Entropy(self.a)+Entropy(self.b)-Entropy(self.m)
205
206    def JaccardInteraction(self):
207        c = Entropy(self.m)
208        if c > 0:
209            return (Entropy(self.a)+Entropy(self.b)-c)/c
210        else:
211            return 0
212
213    def Divergence(self,x,y):
214        ptrue = self.pm[x,y]
215        pkirkwood = self.pa[x]*self.pb[y]
216        if ptrue > 1e-6:
217            div = numpy.log(ptrue/pkirkwood)
218        else:
219            div = 0.0
220        return (ptrue,pkirkwood,_log2*div)
221
222    def ChiSquareP(self):
223        E = numpy.outer(self.pa, self.pb) * self.total
224        return statc.chisqprob(numpy.sum((E-self.m)**2 / E.clip(min=0.000001)), (len(self.pa)-1)*(len(self.pb)-1))
225#        return numpy.sum((E-self.m)**2 / E)
226     
227    def Bootstrap(self,N,limit):
228        # prepare lookup
229        hits = 0
230        nlimit = limit*_log2e
231        f = Flatten(self.m)
232        p = Probabilities(f)
233        LUT = numpy.zeros((self.total,),numpy.int)
234        c = 0
235        for i in xrange(len(f)):
236            for j in xrange(f[i]):
237                LUT[c] = i
238                c += 1
239        assert(c == self.total)
240        for i in xrange(N):
241            nt = numpy.zeros((len(f),),numpy.float)
242            for j in xrange(c):
243                nt[ LUT[randint(0,c-1)] ] += 1
244            q = Probabilities(nt)
245            loss = 0.0
246            for j in xrange(len(f)):
247                if q[j] > 1e-6 and p[j] > 0.0:
248                    loss += q[j]*numpy.log(q[j]/p[j])
249                #loss += p[j]*numpy.log(max(p[j],1e-5)/max(q[j],1e-6))
250            if loss >= nlimit:
251                hits += 1
252        return float(hits)/N
253
254    def Name(self,x,y):
255        s = "%s=%s,%s=%s"%(self.names[0],self.values[0][x],self.names[1],self.values[1][y])
256        return s
257
258    def __init__(self, m, names, values):
259        self.names  = names
260        self.values = values
261        m = numpy.array(m,numpy.float)
262        self.m = m
263
264        self.a = numpy.sum(self.m,axis=1)
265        self.b = numpy.sum(self.m,axis=0)
266        self.total = numpy.sum(self.a)
267
268        self.pa = Probabilities(self.a)
269        self.pb = Probabilities(self.b)
270        self.pm = Probabilities(self.m)
271        dof = 0
272        (ni,nj) = numpy.shape(self.m)
273        for ii in xrange(ni):
274            for jj in xrange(nj):
275                if self.m[ii,jj] > 0:
276                    dof += 1
277        self.dof = dof-1 # degrees of freedom is the number of fields with non-zero counts
278        return
279
280def get3Int(t,a,b,c,wid=None,prior=0):
281    ni = len(a.values)
282    nj = len(b.values)
283    nk = len(c.values)
284    M = prior*numpy.ones((ni,nj,nk),numpy.float)
285    if wid == None: # no weighting
286        for ex in t:
287            if not (ex[a].isSpecial() or ex[b].isSpecial() or ex[c].isSpecial()):
288                M[int(ex[a]),int(ex[b]),int(ex[c])] += 1
289    else:
290        for ex in t:
291            if not (ex[a].isSpecial() or ex[b].isSpecial() or ex[c].isSpecial()):
292                M[int(ex[a]),int(ex[b]),int(ex[c])] += ex.getmeta(wid)
293    N = [a.name,b.name,c.name]
294    V = [[a.values[k] for k in range(ni)],[b.values[k] for k in range(nj)],[c.values[k] for k in range(nk)]]
295    c = ContingencyTable3(M,N,V)
296    return c
297
298def get2Int(t,a,b,wid=None,prior=0):
299    ni = len(a.values)
300    nj = len(b.values)
301    M = prior*numpy.ones((ni,nj),numpy.float)
302    if wid == None: # no weighting
303        for ex in t:
304            if not (ex[a].isSpecial() or ex[b].isSpecial()):
305                M[int(ex[a]),int(ex[b])] += 1
306    else:
307        for ex in t:
308            if not (ex[a].isSpecial() or ex[b].isSpecial()):
309                M[int(ex[a]),int(ex[b])] += ex.getmeta(wid)
310    N = [a.name,b.name]
311    V = [[a.values[k] for k in range(ni)],[b.values[k] for k in range(nj)]]
312    c = ContingencyTable2(M,N,V)
313    return c
314
315def getPvalue(lim,table):
316    # import statistics
317    # return 1-statisticsc.chi_squared(table.dof, 2.0*lim*table.total*_log2e)
318    import statc
319#    print 2.0*lim*table.total*_log2e, table.dof
320    if 2.0*lim*table.total*_log2e <= 0.0:
321        return 1.0 - statc.chisqprob(1.e-10, table.dof)
322    return 1.0 - statc.chisqprob(2.0*lim*table.total*_log2e, table.dof)
323
324def getPvalueDOF(lim,table,dof):
325    # import statisticsc
326    # return 1-statisticsc.chi_squared(dof,2.0*lim*table.total*_log2e)
327    import statc
328    return 1.0 - statc.chisqprob(2.0*lim*table.total*_log2e, dof)
Note: See TracBrowser for help on using the repository browser.