User Tools

Site Tools


pca_class

Visual Analytics Course

PCA (Principal Components Analysis)

Here is a piece of python code that implements the PCA algorithm using a data matrix . You would typically use it as:

data = matrix(...) # some matrix you built or read from disk
 
pca = PCA(data)
pca.computePrincipalComponents()
pca.computeDataProjection() # projects data onto a 2D space

PCA.py
from numpy import *
 
class PCA():
    '''
    a number of functions computing PCA
    and related objects
    vectors a line vectors
    '''
 
 
    def __init__(self, dataMatrix):
        '''
        Constructor
        The method expects a matrix containing all tabular data
        '''
        self.dataMatrix = array(dataMatrix)
        # uno is a 1d vector filled with 1's
        self.uno = array([1] * self.dataMatrix.shape[0])
        '''
            weights is a diagonal matrix with data weights for each data element
            initialized with equal weights
        '''
        self.weights = identity(self.dataMatrix.shape[0], float) * 1.0 / self.dataMatrix.shape[0]
        self.barycenter = None
        self.covarianceMatrix = None
        self.normalizedDataMatrix = None    
 
    def __computeBarycenter__(self):
        self.barycenter = dot(transpose(self.dataMatrix), dot(self.weights, self.uno))
        return True
 
    def __dotVector__(self, v1, v2):
        '''
        v1 and v2 are 1d vectors, the array class dot function would compute
        their scalar product
        array class does not know how to handle the case
        where one wants to multiply a column vector with a row vector
        to output a matrix
        this is what this function does
        '''
 
        m = []
        for i in range(len(v1)):
            for j in range(len(v2)):
                m.append(v1[i] * v2[j])
        m = array(m)
        return m.reshape(len(v1), len(v2))
 
    def __computeCovarianceMatrix__(self):
        '''
        Data normalization
        '''
        if self.barycenter == None:
            self.__computeBarycenter__()
        Y = self.dataMatrix - self.__dotVector__(self.uno, self.barycenter)
        V = dot(dot(transpose(Y), self.weights), Y)
        self.centeredDataMatrix = Y
        d = zeros([self.dataMatrix.shape[1], self.dataMatrix.shape[1]])
        for i in range(self.dataMatrix.shape[1]):
            d[i, i] = 1.0 / std(self.dataMatrix[:, i])
        self.normalizedDataMatrix = dot(Y, d)
        self.covarianceMatrix = dot(dot(d, V), d)
        return True
 
    def computePrincipalComponents(self):
        '''
        Computes *all* eigenValues and eigenVectors of
        the D^T D matrix
        '''
        if self.covarianceMatrix == None:
            self.__computeCovarianceMatrix__()
        self.eigValues, self.eigVectors = linalg.eig(self.covarianceMatrix)
        return True
 
    def computeDataProjection(self):
        '''
        Computes a 2D projection of the data usgin the eigenvectors
        associated with the two largest eigenValues
 
        Assumes method computeCovarianceMatrix has been previously called
        '''
        m = zeros([self.dataMatrix.shape[1], self.dataMatrix.shape[1]])
        for i in range(self.dataMatrix.shape[1]):
            m[i, i] = 1.0 / var(self.dataMatrix[:, i])
        return [dot(dot(self.normalizedDataMatrix, m), self.eigVectors[:, 0]), dot(dot(self.normalizedDataMatrix, m), self.eigVectors[:, 1])]
 
    def cumulativeVariance(self):
        '''
        Computes all ratios
        (lambda_1 + ... + lambda_i)/(lambda_1 + ... + lambda_N)
 
        Looking at the first two values is relevant when projecting in 2D
        '''
        s = 0.0
        for x in self.eigValues:
            s += abs(x)
        v = []
        for x in self.eigValues:
            v.append(abs(x) / s)
        return v
 
    def __project__(self, point, v1 = None, v2 = None):
        '''
        Projects original data points onto the 2D space
        defined by the two given vectors
        When no vector are given, the method uses the
        eigenvectors associated with the two largest eigen values
        '''
        if v1 == None:
            v1 = self.eigVectors[:, 0]
        if v2 == None:
            v2 = self.eigVectors[:, 1]
        return [dot(p, v1), dot(p, v2)]
/net/html/perso/melancon/Visual_Analytics_Course/data/pages/pca_class.txt · Last modified: 2012/10/25 14:57 by melancon