source: orange/Orange/orng/orngCA.py @ 9671:a7b056375472

Revision 9671:a7b056375472, 8.6 KB checked in by anze <anze.staric@…>, 2 years ago (diff)

Moved orange to Orange (part 2)

Line 
1"""
2Correspondence analysis is a descriptive/exploratory technique designed to analyze simple two-way and
3multi-way tables containing some measure of correspondence between the rows and columns.
4"""
5
6from numpy import *
7import numpy.linalg
8import orange
9import operator
10
11def input(filename):
12    """ Loads contingency matrix from the file. """
13    table = [[],[]]
14    try:
15        file = open(filename)
16        try:
17            table = file.readlines()
18        finally:
19            file.close()
20        table = [map(int, el.split()) for el in table]
21    except IOError:
22        pass
23    return table
24   
25class CA(object):
26    """Main class for computation of correspondance analysis"""
27    def __init__(self, contingencyTable, labelR = [], labelC = []):
28        """ @contingencyTable   instance of "list of lists"
29        """     
30        #calculating correspondance analysis from the data matrix
31        #algorithm described in the book (put reference) is used
32
33        self.labelR = labelR
34        self.labelC = labelC
35
36        self.__dataMatrix = matrix(contingencyTable)
37        sumElem = sum(sum(array(self.__dataMatrix))) * 1.
38       
39        #corrMatrix is a matrix of relative frequencies of elements in data matrix
40        self.__corr = self.__dataMatrix / sumElem       
41        self.__colSums = sum(self.__corr, 0)
42        self.__rowSums = sum(self.__corr, 1)
43       
44        self.__colProfiles =  matrix(diag((1. / array(self.__colSums))[0])) * transpose(self.__corr)
45        self.__rowProfiles = matrix(diag((1. / array(self.__rowSums).reshape(1,-1))[0])) * self.__corr
46   
47        self.__a, self.__d, self.__b = self.__calculateSVD();   
48       
49        self.__f = diag((1. / self.__rowSums).reshape(1,-1).tolist()[0]) * self.__a * self.__d
50        self.__g = diag((1. / self.__colSums).tolist()[0]) * self.__b * transpose(self.__d)
51       
52    def __calculateSVD(self):
53        """
54            Computes generalized SVD...
55           
56            This function is used to calculate decomposition A = N * D_mi * M' , where N' * diag(rowSums) * N = I and
57            M' * diag(colSums) * M = I. This decomposition is calculated in 4 steps:
58            i) B = diag(rowSums)^1/2 * A * diag(colSums)^1/2
59            ii) find the ordinary SVD of B: B = U * D * V'
60            iii) N = diag(rowSums)^-1/2 * U
61                M = diag(colSums)^-1/2 * V
62                D_mi = D
63            iv) A= N * D_mi * M'
64           
65            returns (N, D_mi, M)           
66        """
67       
68        a = self.__corr - self.__rowSums * self.__colSums
69        b = diag(sqrt((1. / self.__rowSums).reshape(1,-1).tolist()[0])) * a * diag(sqrt((1. / self.__colSums).tolist()[0]))
70        u, d, v = numpy.linalg.svd(b, 0)
71        N = diag(sqrt(self.__rowSums.reshape(1, -1).tolist()[0])) * u
72        M = diag(sqrt(self.__colSums.tolist()[0])) * transpose(v)
73        d = diag(d.tolist())
74       
75        return (N, d, M)       
76       
77    def getMatrix(self):
78        """
79            Returns array object that is representation of contingency table.
80        """
81        return self.__dataMatrix   
82       
83##    def getCT(self): return self.__ct
84       
85    dataMatrix = property(getMatrix)
86
87    def getA(self): 
88        """
89            columns of A defines the principal axes of the column clouds
90        """
91        return self.__a
92    def getD(self): 
93        """
94            elements on diagonal of D are singular values
95        """
96        return self.__d
97    def getB(self): 
98        """
99            columns of B defines the principal axes of the row clouds
100        """       
101        return self.__b
102    def getF(self): 
103        """
104            coordinates of the row profiles with respect to principal axes B
105        """
106        return self.__f
107    def getG(self): 
108        """
109            coordinates of the column profiles with respect to principal axes A
110        """
111        return self.__g       
112    def getPrincipalRowProfilesCoordinates(self, dim = (0, 1)):
113       """Returns principal co-ordinates of row profiles with respect
114       to principal axis B.
115       Dim defines which principal axes should be taken into account.
116       """
117       if len(dim) == 0:
118           raise Exception("Dim tuple cannot be of lenght zero")
119       return array(take(self.__f, dim, 1))
120    def getPrincipalColProfilesCoordinates(self, dim = (0, 1)): 
121       """Returns principal co-ordinates of column profiles with respect
122       to principal axes A.
123       Dim defines which principal axes should be taken into account.
124       """   
125       if len(dim) == 0:
126           raise Exception("Dim tuple cannot be of lenght zero")     
127       return array(take(self.__g, dim, 1))
128    def getStandardRowCoordinates(self):
129        dinv = where(self.__d != 0, 1. / self.__d, 0)
130        return matrixmultiply(self.__f, transpose(dinv))
131    def getStandardColCoordinates(self):
132        dinv = where(self.__d != 0, 1. / self.__d, 0)
133        return matrixmultiply(self.__g, dinv)
134    def DecompositionOfInertia(self, axis = 0):
135        """
136            axis = 0 decomposition across rows
137            axis = 1 decomposition across columns
138           
139            Columns of this matrix represents contribution of the rows or columns to the inertia of axis.
140        """
141        if axis == 0:
142            return multiply(self.__rowSums, multiply(self.__f, self.__f))
143        else:
144            return multiply(transpose(self.__colSums), multiply(self.__g, self.__g))
145    def InertiaOfAxis(self, percentage = 0):
146        inertias = array(sum(self.DecompositionOfInertia(), 0).tolist()[0])
147        if percentage:
148            return inertias / sum(inertias) * 100
149        else:
150            return inertias
151    def ContributionOfPointsToAxis(self, rowColumn = 0, axis = 0, percentage = 0):
152        contribution = array(transpose(self.DecompositionOfInertia(rowColumn)[:,axis]).tolist()[0])
153        if percentage:
154            return contribution / sum(contribution) * 100
155        else:
156            return contribution
157    def PointsWithMostInertia(self, rowColumn = 0, axis = (0, 1)):
158        contribution = self.ContributionOfPointsToAxis(rowColumn = rowColumn, axis = axis[0], percentage = 0) + \
159                        self.ContributionOfPointsToAxis(rowColumn = rowColumn, axis = axis[1], percentage = 0)
160        tmp = zip(range(len(contribution)), contribution)
161
162        tmp.sort(lambda x, y: cmp(x[1], y[1]))
163
164        a = [i for (i, v) in tmp]
165        a.reverse()
166        return a
167   
168    def PlotScreeDiagram(self):
169        # todo: legend, axis, etc
170        import pylab
171        pylab.plot(range(1, min(self.__dataMatrix.shape) + 1), self.InertiaOfAxis(1))
172        pylab.axis([0, min(self.__dataMatrix.shape) + 1, 0, 100])
173        pylab.show()
174       
175    def Biplot(self, dim = (0, 1)):
176        import pylab
177        if len(dim) != 2:
178           raise Exception("Dim tuple must be of length two")
179        pylab.plot(self.getPrincipalRowProfilesCoordinates()[:, dim[0]], self.getPrincipalRowProfilesCoordinates()[:, dim[1]], 'ro',
180            self.getPrincipalColProfilesCoordinates()[:, dim[0]], self.getPrincipalColProfilesCoordinates()[:, dim[1]], 'bs')
181        if self.labelR:
182            for i, x, y in zip(range(len(self.getPrincipalRowProfilesCoordinates()[:, dim[0]])), \
183                                    self.getPrincipalRowProfilesCoordinates()[:, dim[0]], \
184                                    self.getPrincipalRowProfilesCoordinates()[:, dim[1]]):
185                pylab.text(x, y, self.labelR[i], horizontalalignment='center')
186        if self.labelC:
187            for i, x, y in zip(range(len(self.getPrincipalColProfilesCoordinates()[:, dim[0]])), \
188                                    self.getPrincipalColProfilesCoordinates()[:, dim[0]], \
189                                    self.getPrincipalColProfilesCoordinates()[:, dim[1]]):
190                pylab.text(x, y, self.labelC[i], horizontalalignment='center')               
191        pylab.grid()
192        pylab.show()               
193   
194    A = property(getA)
195    B = property(getB)
196    D = property(getD)
197    F = property(getF)
198    G = property(getG)
199   
200if __name__ == '__main__':
201##    a = random.random_integers(0, 100, 100).reshape(10,-1)
202##    c = CA(a)
203##    c.Biplot()
204
205##    data = matrix([[72,    39,    26,    23 ,    4],
206##    [95,    58,    66,    84,    41],
207##    [80,    73,    83,     4 ,   96],
208##    [79,    93,    35,    73,    63]])
209##
210##    data = [[9, 11, 4],
211##                [ 3,          5,          3],
212##                [     11,          6,          3],
213##                [24,         73,         48]]
214
215   
216    data = input('doc\\datasets\\smokers.tab')
217    c = CA(data, ['Senior Managers', 'Junior Managers', 'Senior Employees', 'Junior Employees', 'Secretaries'], 
218        ['None', 'Light', 'Medium', 'Heavy'])
219#    c.PlotScreeDiagram()
Note: See TracBrowser for help on using the repository browser.