Ignore:
Timestamp:
03/26/12 16:37:25 (2 years ago)
Author:
anze <anze.staric@…>
Branch:
default
Message:

PCA now returns the same results as prcomp(X, scale=TRUE) in R.
(Fixed normalization in linalg numpy operations and output of st. dev)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • Orange/projection/linear.py

    r10644 r10645  
    163163 
    164164        return Projector(input_domain=domain, 
    165                          mean=Xm, 
    166                          stdev=stdev, 
     165                         center=Xm, 
     166                         scale=stdev, 
    167167                         standardize=False, 
    168168                         projection=U) 
     
    933933        if method == DR_PCA: 
    934934            pca = Pca(standardize=False, max_components=ncomps, 
    935                       use_generalized_eigenvectors=False) 
     935                      use_generalized_eigenvectors=False, ddof=0) 
    936936            domain = data.Domain([feature.Continuous("g%d" % i) for i 
    937937                                  in xrange(len(data_matrix))], False) 
     
    940940        elif method == DR_SPCA and self.graph.data_has_class: 
    941941            pca = Spca(standardize=False, max_components=ncomps, 
    942                        use_generalized_eigenvectors=self.use_generalized_eigenvectors) 
     942                       use_generalized_eigenvectors=self.use_generalized_eigenvectors, ddof=0) 
    943943            domain = data.Domain([feature.Continuous("g%d" % i) for i 
    944944                                  in xrange(len(data_matrix))], feature.Continuous("c")) 
     
    12091209    #: Array containing means of each variable in the data set that was used 
    12101210    #: to construct the projection. 
    1211     mean = numpy.array(()) 
     1211    center = numpy.array(()) 
    12121212 
    12131213    #: An array containing standard deviations of each variable in the data 
    12141214    #: set that was used to construct the projection. 
    1215     stdev = numpy.array(()) 
     1215    scale = numpy.array(()) 
    12161216 
    12171217    #: True, if standardization was used when constructing the projection. If 
     
    12381238    def _project_single(self, example, return_what, new_feature, feature_idx): 
    12391239        ex = Orange.data.Table([example]).to_numpy("a")[0] 
    1240         ex -= self.mean 
     1240        ex -= self.center 
    12411241        if self.standardize: 
    1242             ex /= self.stdev 
     1242            ex /= self.scale 
    12431243        return new_feature(numpy.dot(self.projection[feature_idx, :], ex.T)[0]) 
    12441244 
     
    12581258 
    12591259        X, = dataset.to_numpy_MA("a") 
    1260         Xm, U = self.mean, self.projection 
     1260        Xm, U = self.center, self.projection 
    12611261        n, m = X.shape 
    12621262 
     
    12671267 
    12681268        if self.standardize: 
    1269             Xd /= self.stdev 
     1269            Xd /= self.scale 
    12701270 
    12711271        self.A = numpy.dot(Xd, U.T) 
     
    12951295    """ 
    12961296 
     1297    #: Delta degrees of freedom used for numpy operations. 
     1298    #: 1 means normalization with (N-1) in cov and std operations 
     1299    ddof = 1 
     1300 
    12971301    def __new__(cls, dataset=None, **kwds): 
    12981302        optimizer = object.__new__(cls) 
     
    13401344 
    13411345        return PcaProjector(input_domain=dataset.domain, 
    1342                             mean=mean, 
    1343                             stdev=stdev, 
     1346                            center=mean, 
     1347                            scale=stdev, 
    13441348                            standardize=self.standardize, 
    1345                             eigen_vectors=components, 
    13461349                            projection=components, 
    1347                             eigen_values=variances, 
     1350                            variances=variances, 
    13481351                            variance_sum=variance_sum) 
    13491352 
     
    13601363 
    13611364        #take care of the constant features 
    1362         stdev = numpy.std(Xd, axis=0) 
     1365        stdev = numpy.std(Xd, axis=0, ddof=self.ddof) 
    13631366        stdev[stdev == 0] = 1. # to prevent division by zero 
    13641367        relevant_features = stdev != 0 
     
    13711374        n, m = Xd.shape 
    13721375        if n < m: 
    1373             C = numpy.dot(Xg.T, Xd.T) 
     1376            C = numpy.dot(Xg.T, Xd.T) / (m - self.ddof) 
    13741377            V, D, T = numpy.linalg.svd(C) 
    1375             U = numpy.dot(V.T, Xd) / numpy.sqrt(D.reshape(-1, 1)) 
     1378            U = numpy.dot(V.T, Xd) / numpy.sqrt(D.reshape(-1, 1) * (m - self.ddof)) 
    13761379        else: 
    1377             C = numpy.dot(Xg, Xd) 
     1380            C = numpy.dot(Xg, Xd) / (n - self.ddof) 
    13781381            U, D, T = numpy.linalg.svd(C) 
    13791382            U = U.T  # eigenvectors are now in rows 
     
    14041407 
    14051408 
    1406 class Spca(Pca): 
     1409class Spca(PCA): 
    14071410    def _perform_pca(self, dataset, Xd, Xg): 
    14081411        # define the Laplacian matrix 
     
    14181421@deprecated_members({"pc_domain": "output_domain"}) 
    14191422class PcaProjector(Projector): 
    1420     #: Synonymous for :obj:`~Orange.projection.linear.Projector.projection`. 
    1421     eigen_vectors = numpy.array(()).reshape(0, 0) 
    1422  
    1423     #: Array containing standard deviations of principal components. 
    1424     eigen_values = numpy.array(()) 
     1423    #: Array containing variances of principal components. 
     1424    variances = numpy.array(()) 
    14251425 
    14261426    #: Sum of all variances in the data set that was used to construct the PCA space. 
     
    14301430        ncomponents = 10 
    14311431        s = self.variance_sum 
    1432         cs = numpy.cumsum(self.eigen_values) / s 
     1432        cs = numpy.cumsum(self.variances) / s 
     1433        stdev = numpy.sqrt(self.variances) 
    14331434        return "\n".join([ 
    14341435            "PCA SUMMARY", 
     
    14381439                     ["%10s" % a.name for a in self.output_domain.attributes]), 
    14391440            " ".join(["Std. deviation"] + 
    1440                      ["%10.3f" % a for a in self.eigen_values]), 
     1441                     ["%10.3f" % a for a in stdev]), 
    14411442            " ".join(["Proportion Var"] + 
    1442                      ["%10.3f" % a for a in self.eigen_values / s * 100]), 
     1443                     ["%10.3f" % a for a in self.variances / s * 100]), 
    14431444            " ".join(["Cumulative Var"] + 
    14441445                     ["%10.3f" % a for a in cs * 100]), 
     
    14541455                     ["%10s" % self.output_domain.attributes[-1].name]), 
    14551456            " ".join(["Std. deviation"] + 
    1456                      ["%10.3f" % a for a in self.eigen_values[:ncomponents]] + 
     1457                     ["%10.3f" % a for a in stdev[:ncomponents]] + 
    14571458                     ["%10s" % ""] + 
    1458                      ["%10.3f" % self.eigen_values[-1]]), 
     1459                     ["%10.3f" % stdev[-1]]), 
    14591460            " ".join(["Proportion Var"] + 
    1460                      ["%10.3f" % a for a in self.eigen_values[:ncomponents] / s * 100] + 
     1461                     ["%10.3f" % a for a in self.variances[:ncomponents] / s * 100] + 
    14611462                     ["%10s" % ""] + 
    1462                      ["%10.3f" % (self.eigen_values[-1] / s * 100)]), 
     1463                     ["%10.3f" % (self.variances[-1] / s * 100)]), 
    14631464            " ".join(["Cumulative Var"] + 
    14641465                     ["%10.3f" % a for a in cs[:ncomponents] * 100] + 
     
    14821483 
    14831484        s = self.variance_sum 
    1484         vc = self.eigen_values / s 
    1485         cs = numpy.cumsum(self.eigen_values) / s 
     1485        vc = self.variances / s 
     1486        cs = numpy.cumsum(self.variances) / s 
    14861487 
    14871488        fig = plt.figure() 
    14881489        ax = fig.add_subplot(111) 
    14891490 
    1490         x_axis = range(len(self.eigen_values)) 
     1491        x_axis = range(len(self.variances)) 
    14911492        plt.grid(True) 
    14921493 
     
    15011502        ax.legend(loc=0) 
    15021503 
    1503         ax.axis([-0.5, len(self.eigen_values) - 0.5, 0, 1]) 
     1504        ax.axis([-0.5, len(self.variances) - 0.5, 0, 1]) 
    15041505 
    15051506        if filename: 
     
    15261527            raise ValueError, 'Two components are needed for biplot' 
    15271528 
    1528         if not (0 <= min(components) <= max(components) < len(self.eigen_values)): 
     1529        if not (0 <= min(components) <= max(components) < len(self.variances)): 
    15291530            raise ValueError, 'Invalid components' 
    15301531 
     
    15321533        Y = self.A[:, components[1]] 
    15331534 
    1534         vectorsX = self.eigen_vectors[:, components[0]] 
    1535         vectorsY = self.eigen_vectors[:, components[1]] 
    1536  
    1537         max_load_value = self.eigen_vectors.max() * 1.5 
     1535        vectorsX = self.variances[:, components[0]] 
     1536        vectorsY = self.variances[:, components[1]] 
     1537 
     1538        max_load_value = self.variances.max() * 1.5 
    15381539 
    15391540        fig = plt.figure() 
    15401541        ax1 = fig.add_subplot(111) 
    15411542        ax1.set_title(title + "\n") 
    1542         ax1.set_xlabel("PC%s (%d%%)" % (components[0], self.eigen_values[components[0]] / self.variance_sum * 100)) 
    1543         ax1.set_ylabel("PC%s (%d%%)" % (components[1], self.eigen_values[components[1]] / self.variance_sum * 100)) 
     1543        ax1.set_xlabel("PC%s (%d%%)" % (components[0], self.variances[components[0]] / self.variance_sum * 100)) 
     1544        ax1.set_ylabel("PC%s (%d%%)" % (components[1], self.variances[components[1]] / self.variance_sum * 100)) 
    15441545        ax1.xaxis.set_label_position('bottom') 
    15451546        ax1.xaxis.set_ticks_position('bottom') 
     
    16681669        return FdaProjector(input_domain=dataset.domain, 
    16691670                            output_domain=out_domain, 
    1670                             mean=Xm, 
    1671                             stdev=stdev, 
     1671                            center=Xm, 
     1672                            scale=stdev, 
    16721673                            standardize=True, 
    16731674                            eigen_vectors=U, 
Note: See TracChangeset for help on using the changeset viewer.