source: orange/Orange/projection/pca.py @ 9927:d6ca7b346864

Revision 9927:d6ca7b346864, 12.9 KB checked in by markotoplak, 2 years ago (diff)

data.variable -> feature.

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, dataset):
148        if type(dataset) != Orange.data.Table:
149            dataset = Orange.data.Table([dataset])
150
151        X = dataset.to_numpy_MA("a")[0]
152        Xm, U = self.mean, self.eigen_vectors
153        n, m = X.shape
154       
155        if m != len(self.stdev):
156            raise orange.KernelException, "Invalid number of features"
157       
158        Xd = X - Xm
159       
160        if self.standardize:
161            Xd /= self.stdev
162       
163        self.A = np.ma.dot(Xd, U.T)
164       
165        return Orange.data.Table(self.pc_domain, self.A.tolist())
166   
167    def __str__(self):
168        ncomponents = 10
169        s = self.variance_sum
170        cs = np.cumsum(self.eigen_values) / s
171        return "\n".join([
172        "PCA SUMMARY",
173        "",
174        "Std. deviation of components:",
175        " ".join(["              "] +
176                 ["%10s" % a.name for a in self.pc_domain.attributes]),
177        " ".join(["Std. deviation"] +
178                 ["%10.3f" % a for a in self.eigen_values]),
179        " ".join(["Proportion Var"] + 
180                 ["%10.3f" % a for a in  self.eigen_values / s * 100]),
181        " ".join(["Cumulative Var"] +
182                 ["%10.3f" % a for a in cs * 100]),
183        "",
184        #"Loadings:",
185        #" ".join(["%10s"%""] + ["%10s" % a.name for a in self.pc_domain]),
186        #"\n".join([
187        #    " ".join([a.name] + ["%10.3f" % b for b in self.eigen_vectors.T[i]])
188        #          for i, a in enumerate(self.input_domain.attributes)
189        #          ])
190        ]) if len(self.pc_domain) <= ncomponents else \
191        "\n".join([
192        "PCA SUMMARY",
193        "",
194        "Std. deviation of components:",
195        " ".join(["              "] +
196                 ["%10s" % a.name for a in self.pc_domain.attributes[:ncomponents]] +
197                 ["%10s" % "..."] +
198                 ["%10s" % self.pc_domain.attributes[-1].name]),
199        " ".join(["Std. deviation"] +
200                 ["%10.3f" % a for a in self.eigen_values[:ncomponents]] + 
201                 ["%10s" % ""] +
202                 ["%10.3f" % self.eigen_values[-1]]),
203        " ".join(["Proportion Var"] + 
204                 ["%10.3f" % a for a in self.eigen_values[:ncomponents] / s * 100] + 
205                 ["%10s" % ""] +
206                 ["%10.3f" % (self.eigen_values[-1] / s * 100)]),
207        " ".join(["Cumulative Var"] +
208                 ["%10.3f" % a for a in cs[:ncomponents] * 100] + 
209                 ["%10s" % ""] +
210                 ["%10.3f" % (cs[-1] * 100)]),
211        "",
212        #"Loadings:",
213        #" ".join(["%16s" % ""] +
214        #         ["%8s" % a.name for a in self.pc_domain.attributes[:ncomponents]] +
215        #         ["%8s" % "..."] +
216        #         ["%8s" % self.pc_domain.attributes[-1].name]),
217        #"\n".join([
218        #    " ".join(["%16.16s" %a.name] +
219        #             ["%8.3f" % b for b in self.eigen_vectors.T[i, :ncomponents]] +
220        #             ["%8s" % ""] +
221        #             ["%8.3f" % self.eigen_vectors.T[i, -1]])
222        #          for i, a in enumerate(self.input_domain.attributes)
223        #          ])
224        ])
225
226       
227   
228    ################ Plotting functions ###################
229   
230    def scree_plot(self, filename = None, title = 'Scree plot'):
231        """
232        Draw a scree plot of principal components
233       
234        :param filename: Name of the file to which the plot will be saved. \
235        If None, plot will be displayed instead.
236        :type filename: str
237        :param title: Plot title
238        :type title: str
239        """
240        import pylab as plt
241       
242        s = self.variance_sum
243        vc = self.eigen_values / s
244        cs = np.cumsum(self.eigen_values) / s
245       
246        fig = plt.figure()
247        ax = fig.add_subplot(111)
248       
249        x_axis = range(len(self.eigen_values))
250        x_labels = ["PC%d" % (i + 1, ) for i in x_axis]
251       
252        ax.set_xticks(x_axis)
253        ax.set_xticklabels(x_labels)
254        plt.setp(ax.get_xticklabels(), "rotation", 90)
255        plt.grid(True)
256       
257        ax.set_xlabel('Principal components')
258        ax.set_ylabel('Proportion of Variance')
259        ax.set_title(title + "\n")
260        ax.plot(x_axis, vc, color="red")
261        ax.scatter(x_axis, vc, color="red", label="Variance")
262       
263        ax.plot(x_axis, cs, color="orange")
264        ax.scatter(x_axis, cs, color="orange", label="Cumulative Variance")
265        ax.legend(loc=0)
266       
267        ax.axis([-0.5, len(self.eigen_values) - 0.5, 0, 1])
268       
269        if filename:
270            plt.savefig(filename)
271        else:
272            plt.show()
273           
274    def biplot(self, filename = None, components = [0,1], title = 'Biplot'):
275        """
276        Draw biplot for PCA. Actual projection must be performed via pca(data)
277        before bipot can be used.
278       
279        :param filename: Name of the file to which the plot will be saved. \
280        If None, plot will be displayed instead.
281        :type plot: str
282        :param components: List of two components to plot.
283        :type components: list
284        :param title: Plot title
285        :type title: str
286        """
287        import pylab as plt
288       
289        if len(components) < 2:
290            raise orange.KernelException, 'Two components are needed for biplot'
291       
292        if not (0 <= min(components) <= max(components) < len(self.eigen_values)):
293            raise orange.KernelException, 'Invalid components'
294       
295        X = self.A[:,components[0]]
296        Y = self.A[:,components[1]]
297       
298        vectorsX = self.eigen_vectors[:,components[0]]
299        vectorsY = self.eigen_vectors[:,components[1]]
300       
301       
302        #TO DO -> pc.biplot (maybe)
303        #trDataMatrix = dataMatrix / lam
304        #trLoadings = loadings * lam
305       
306        #max_data_value = numpy.max(abs(trDataMatrix)) * 1.05
307        max_load_value = self.eigen_vectors.max() * 1.5
308       
309        #plt.clf()
310        fig = plt.figure()
311        ax1 = fig.add_subplot(111)
312        ax1.set_title(title + "\n")
313        ax1.set_xlabel("PC%s (%d%%)" % (components[0], self.eigen_values[components[0]] / self.variance_sum * 100))
314        ax1.set_ylabel("PC%s (%d%%)" % (components[1], self.eigen_values[components[1]] / self.variance_sum * 100))
315        ax1.xaxis.set_label_position('bottom')
316        ax1.xaxis.set_ticks_position('bottom')
317        ax1.yaxis.set_label_position('left')
318        ax1.yaxis.set_ticks_position('left')       
319       
320        #if self._classArray == None:
321        #trDataMatrix = transpose(trDataMatrix)
322        ax1.plot(X, Y, Colors[0])
323        #else:
324            #suboptimal
325        #    classValues = []
326        #    for classValue in self._classArray:
327        #        if classValue not in classValues:
328        #            classValues.append(classValue)
329        #    for i in range(len(classValues)):
330        #        choice = numpy.array([classValues[i] == cv for cv in self._classArray])
331        #        partialDataMatrix = transpose(trDataMatrix[choice])
332        #        ax1.plot(partialDataMatrix[0], partialDataMatrix[1],
333        #                 Colors[i % len(Colors)], label = str(classValues[i]))
334        #    ax1.legend()
335       
336        #ax1.set_xlim(-max_data_value, max_data_value)
337        #ax1.set_ylim(-max_data_value, max_data_value)
338       
339        #eliminate double axis on right
340        ax0 = ax1.twinx()
341        ax0.yaxis.set_visible(False)
342               
343        ax2 = ax0.twiny()
344        ax2.xaxis.set_label_position('top')
345        ax2.xaxis.set_ticks_position('top')
346        ax2.yaxis.set_label_position('right')
347        ax2.yaxis.set_ticks_position('right')
348        for tl in ax2.get_xticklabels():
349            tl.set_color('r')
350        for tl in ax2.get_yticklabels():
351            tl.set_color('r')
352       
353        arrowprops = dict(facecolor = 'red', edgecolor = 'red', width = 1, headwidth = 4)
354
355        for (x, y, a) in zip(vectorsX, vectorsY,self.input_domain.attributes):
356            if max(x, y) < 0.1:
357                continue
358            print x, y, a
359            ax2.annotate('', (x, y), (0, 0), arrowprops = arrowprops)
360            ax2.text(x * 1.1, y * 1.2, a.name, color = 'red')
361           
362        ax2.set_xlim(-max_load_value, max_load_value)
363        ax2.set_ylim(-max_load_value, max_load_value)
364       
365        if filename:
366            plt.savefig(filename)
367        else:
368            plt.show()
369 
Note: See TracBrowser for help on using the repository browser.