source: orange/Orange/projection/pca.py @ 10181:b839b8b6fc00

Revision 10181:b839b8b6fc00, 13.0 KB checked in by anze <anze.staric@…>, 2 years ago (diff)

Added docstring for call method in PcaClassifier.

Line 
1import Orange.data
2import Orange.feature
3import numpy as np
4   
5#color table for biplot
6Colors = ['bo','go','yo','co','mo']
7
8class Pca(object):
9    """
10    Orthogonal transformation of data into a set of uncorrelated variables called
11    principal components. This transformation is defined in such a way that the
12    first variable has as high variance as possible.
13
14    If data instances are provided to the constructor, the learning algorithm
15    is called and the resulting classifier is returned instead of the learner.
16
17    :param standardize: Perform standardization of the dataset
18    :type standardize: boolean
19    :param max_components: Maximum number of retained components
20    :type max_components: int
21    :param variance_covered: Percent of the variance to cover with components
22    :type variance_covered: float
23       
24    :rtype: :class:`Orange.projection.pca.Pca` or
25            :class:`Orange.projection.pca.PcaClassifier`
26    """
27   
28    def __new__(cls, dataset = None, **kwds):
29        learner = object.__new__(cls)
30        learner.__init__(**kwds)
31       
32        if dataset:
33            return learner(dataset)
34        else:
35            return learner
36   
37    def __init__(self, standardize = True,
38                 max_components = 0, variance_covered = 1):
39        self.standardize = standardize
40        self.max_components = max_components
41        self.variance_covered = variance_covered if variance_covered < 1. else 1
42   
43    def __call__(self, dataset):
44        """
45        Perform a pca analysis on a dataset and return a classifier that maps data
46        into principal component subspace.
47        """
48       
49        X = dataset.to_numpy_MA("a")[0]
50        N,M = X.shape
51        Xm = np.mean(X, axis=0)
52        Xd = X - Xm
53       
54        #take care of the constant features
55        stdev = np.std(Xd, axis=0)
56        relevant_features = stdev != 0
57        if self.standardize:
58            stdev[stdev == 0] = 1.
59            Xd /= stdev
60        Xd = Xd[:,relevant_features]
61       
62        #actual pca
63        n,m = Xd.shape
64        if n < m:
65            C = np.ma.dot(Xd, Xd.T)
66            V, D, T = np.linalg.svd(C)
67            U = np.ma.dot(V.T, Xd) / np.sqrt(D.reshape(-1,1))
68        else:
69            C = np.ma.dot(Xd.T, Xd)
70            U, D, T = np.linalg.svd(C)
71       
72        #insert zeros for constant features
73        n, m = U.shape
74        if m != M:
75            U_ = np.zeros((n,M))
76            U_[:,relevant_features] = U
77            U = U_
78       
79        variance_sum = D.sum()
80       
81        #select eigen vectors
82        if self.variance_covered != 1:
83            nfeatures = np.nonzero(np.cumsum(D) / sum(D) >= self.variance_covered)[0][0] + 1
84            U = U[:nfeatures, :]
85            D = D[:nfeatures]
86       
87        if self.max_components > 0:
88            U = U[:self.max_components, :]
89            D = D[:self.max_components]
90       
91        n, m = U.shape
92        pc_domain = Orange.data.Domain([Orange.feature.Continuous("Comp.%d"%(i+1)) for i in range(n)], False)
93       
94        return PcaClassifier(input_domain = dataset.domain,
95                             pc_domain = pc_domain,
96                             mean = Xm,
97                             stdev = stdev,
98                             standardize = self.standardize,
99                             eigen_vectors = U,
100                             eigen_values = D,
101                             variance_sum = variance_sum)
102
103
104class PcaClassifier(object):
105    """
106    .. attribute:: input_domain
107   
108        Domain of the dataset that was used to construct principal component
109        subspace.
110       
111    .. attribute:: pc_domain
112   
113        Domain used in returned datasets. This domain has a Float variable for
114        each principal component and no class variable.
115       
116    .. attribute:: mean
117   
118        Array containing means of each variable in the dataset that was used
119        to construct pca space.
120       
121    .. attribute:: stdev
122   
123        An array containing standard deviations of each variable in the dataset
124        that was used to construct pca space.
125       
126    .. attribute:: standardize
127   
128        True, if standardization was used when constructing the pca space. If set,
129        instances will be standardized before being mapped to the pca space.
130   
131    .. attribute:: eigen_vectors
132   
133        Array containing vectors that are used to map to pca space.
134       
135    .. attribute:: eigen_values
136   
137        Array containing standard deviations of principal components.
138   
139    .. attribute:: variance_sum
140   
141        Sum of all variances in the dataset that was used to construct the pca
142        space.
143    """
144    def __init__(self, **kwds):
145        self.__dict__.update(kwds)
146   
147    def __call__(self, data):
148        """
149        Project data into PCA feature space.
150
151        :param data: input dataset
152
153        :rtype: :obj:`~Orange.data.Table`
154        """
155        if type(dataset) != Orange.data.Table:
156            dataset = Orange.data.Table([dataset])
157
158        X = dataset.to_numpy_MA("a")[0]
159        Xm, U = self.mean, self.eigen_vectors
160        n, m = X.shape
161       
162        if m != len(self.stdev):
163            raise orange.KernelException, "Invalid number of features"
164       
165        Xd = X - Xm
166       
167        if self.standardize:
168            Xd /= self.stdev
169       
170        self.A = np.ma.dot(Xd, U.T)
171       
172        return Orange.data.Table(self.pc_domain, self.A.tolist())
173   
174    def __str__(self):
175        ncomponents = 10
176        s = self.variance_sum
177        cs = np.cumsum(self.eigen_values) / s
178        return "\n".join([
179        "PCA SUMMARY",
180        "",
181        "Std. deviation of components:",
182        " ".join(["              "] +
183                 ["%10s" % a.name for a in self.pc_domain.attributes]),
184        " ".join(["Std. deviation"] +
185                 ["%10.3f" % a for a in self.eigen_values]),
186        " ".join(["Proportion Var"] + 
187                 ["%10.3f" % a for a in  self.eigen_values / s * 100]),
188        " ".join(["Cumulative Var"] +
189                 ["%10.3f" % a for a in cs * 100]),
190        "",
191        #"Loadings:",
192        #" ".join(["%10s"%""] + ["%10s" % a.name for a in self.pc_domain]),
193        #"\n".join([
194        #    " ".join([a.name] + ["%10.3f" % b for b in self.eigen_vectors.T[i]])
195        #          for i, a in enumerate(self.input_domain.attributes)
196        #          ])
197        ]) if len(self.pc_domain) <= ncomponents else \
198        "\n".join([
199        "PCA SUMMARY",
200        "",
201        "Std. deviation of components:",
202        " ".join(["              "] +
203                 ["%10s" % a.name for a in self.pc_domain.attributes[:ncomponents]] +
204                 ["%10s" % "..."] +
205                 ["%10s" % self.pc_domain.attributes[-1].name]),
206        " ".join(["Std. deviation"] +
207                 ["%10.3f" % a for a in self.eigen_values[:ncomponents]] + 
208                 ["%10s" % ""] +
209                 ["%10.3f" % self.eigen_values[-1]]),
210        " ".join(["Proportion Var"] + 
211                 ["%10.3f" % a for a in self.eigen_values[:ncomponents] / s * 100] + 
212                 ["%10s" % ""] +
213                 ["%10.3f" % (self.eigen_values[-1] / s * 100)]),
214        " ".join(["Cumulative Var"] +
215                 ["%10.3f" % a for a in cs[:ncomponents] * 100] + 
216                 ["%10s" % ""] +
217                 ["%10.3f" % (cs[-1] * 100)]),
218        "",
219        #"Loadings:",
220        #" ".join(["%16s" % ""] +
221        #         ["%8s" % a.name for a in self.pc_domain.attributes[:ncomponents]] +
222        #         ["%8s" % "..."] +
223        #         ["%8s" % self.pc_domain.attributes[-1].name]),
224        #"\n".join([
225        #    " ".join(["%16.16s" %a.name] +
226        #             ["%8.3f" % b for b in self.eigen_vectors.T[i, :ncomponents]] +
227        #             ["%8s" % ""] +
228        #             ["%8.3f" % self.eigen_vectors.T[i, -1]])
229        #          for i, a in enumerate(self.input_domain.attributes)
230        #          ])
231        ])
232
233       
234   
235    ################ Plotting functions ###################
236   
237    def scree_plot(self, filename = None, title = 'Scree plot'):
238        """
239        Draw a scree plot of principal components
240       
241        :param filename: Name of the file to which the plot will be saved. \
242        If None, plot will be displayed instead.
243        :type filename: str
244        :param title: Plot title
245        :type title: str
246        """
247        import pylab as plt
248       
249        s = self.variance_sum
250        vc = self.eigen_values / s
251        cs = np.cumsum(self.eigen_values) / s
252       
253        fig = plt.figure()
254        ax = fig.add_subplot(111)
255       
256        x_axis = range(len(self.eigen_values))
257        x_labels = ["PC%d" % (i + 1, ) for i in x_axis]
258       
259        ax.set_xticks(x_axis)
260        ax.set_xticklabels(x_labels)
261        plt.setp(ax.get_xticklabels(), "rotation", 90)
262        plt.grid(True)
263       
264        ax.set_xlabel('Principal components')
265        ax.set_ylabel('Proportion of Variance')
266        ax.set_title(title + "\n")
267        ax.plot(x_axis, vc, color="red")
268        ax.scatter(x_axis, vc, color="red", label="Variance")
269       
270        ax.plot(x_axis, cs, color="orange")
271        ax.scatter(x_axis, cs, color="orange", label="Cumulative Variance")
272        ax.legend(loc=0)
273       
274        ax.axis([-0.5, len(self.eigen_values) - 0.5, 0, 1])
275       
276        if filename:
277            plt.savefig(filename)
278        else:
279            plt.show()
280           
281    def biplot(self, filename = None, components = [0,1], title = 'Biplot'):
282        """
283        Draw biplot for PCA. Actual projection must be performed via pca(data)
284        before bipot can be used.
285       
286        :param filename: Name of the file to which the plot will be saved. \
287        If None, plot will be displayed instead.
288        :type plot: str
289        :param components: List of two components to plot.
290        :type components: list
291        :param title: Plot title
292        :type title: str
293        """
294        import pylab as plt
295       
296        if len(components) < 2:
297            raise orange.KernelException, 'Two components are needed for biplot'
298       
299        if not (0 <= min(components) <= max(components) < len(self.eigen_values)):
300            raise orange.KernelException, 'Invalid components'
301       
302        X = self.A[:,components[0]]
303        Y = self.A[:,components[1]]
304       
305        vectorsX = self.eigen_vectors[:,components[0]]
306        vectorsY = self.eigen_vectors[:,components[1]]
307       
308       
309        #TO DO -> pc.biplot (maybe)
310        #trDataMatrix = dataMatrix / lam
311        #trLoadings = loadings * lam
312       
313        #max_data_value = numpy.max(abs(trDataMatrix)) * 1.05
314        max_load_value = self.eigen_vectors.max() * 1.5
315       
316        #plt.clf()
317        fig = plt.figure()
318        ax1 = fig.add_subplot(111)
319        ax1.set_title(title + "\n")
320        ax1.set_xlabel("PC%s (%d%%)" % (components[0], self.eigen_values[components[0]] / self.variance_sum * 100))
321        ax1.set_ylabel("PC%s (%d%%)" % (components[1], self.eigen_values[components[1]] / self.variance_sum * 100))
322        ax1.xaxis.set_label_position('bottom')
323        ax1.xaxis.set_ticks_position('bottom')
324        ax1.yaxis.set_label_position('left')
325        ax1.yaxis.set_ticks_position('left')       
326       
327        #if self._classArray == None:
328        #trDataMatrix = transpose(trDataMatrix)
329        ax1.plot(X, Y, Colors[0])
330        #else:
331            #suboptimal
332        #    classValues = []
333        #    for classValue in self._classArray:
334        #        if classValue not in classValues:
335        #            classValues.append(classValue)
336        #    for i in range(len(classValues)):
337        #        choice = numpy.array([classValues[i] == cv for cv in self._classArray])
338        #        partialDataMatrix = transpose(trDataMatrix[choice])
339        #        ax1.plot(partialDataMatrix[0], partialDataMatrix[1],
340        #                 Colors[i % len(Colors)], label = str(classValues[i]))
341        #    ax1.legend()
342       
343        #ax1.set_xlim(-max_data_value, max_data_value)
344        #ax1.set_ylim(-max_data_value, max_data_value)
345       
346        #eliminate double axis on right
347        ax0 = ax1.twinx()
348        ax0.yaxis.set_visible(False)
349               
350        ax2 = ax0.twiny()
351        ax2.xaxis.set_label_position('top')
352        ax2.xaxis.set_ticks_position('top')
353        ax2.yaxis.set_label_position('right')
354        ax2.yaxis.set_ticks_position('right')
355        for tl in ax2.get_xticklabels():
356            tl.set_color('r')
357        for tl in ax2.get_yticklabels():
358            tl.set_color('r')
359       
360        arrowprops = dict(facecolor = 'red', edgecolor = 'red', width = 1, headwidth = 4)
361
362        for (x, y, a) in zip(vectorsX, vectorsY,self.input_domain.attributes):
363            if max(x, y) < 0.1:
364                continue
365            print x, y, a
366            ax2.annotate('', (x, y), (0, 0), arrowprops = arrowprops)
367            ax2.text(x * 1.1, y * 1.2, a.name, color = 'red')
368           
369        ax2.set_xlim(-max_load_value, max_load_value)
370        ax2.set_ylim(-max_load_value, max_load_value)
371       
372        if filename:
373            plt.savefig(filename)
374        else:
375            plt.show()
376 
Note: See TracBrowser for help on using the repository browser.