Python 主成分分析PCA

http://www.cnblogs.com/chenbjin/p/4200790.html

　　主成分分析（PCA）是一種基於變數協方差矩陣對數據進行壓縮降維、去噪的有效方法，PCA的思想是將n維特徵映射到k維上（k<n），這k維特徵稱為主元，是舊特徵的線性組合，這些線性組合最大化樣本方差，盡量使新的k個特徵互不相關。

相關知識

介紹一個PCA的教程：A tutorial on Principal Components Analysis ——Lindsay I Smith

1.協方差Covariance

　　變數X和變數Y的協方差公式如下，協方差是描述不同變數之間的相關關係，協方差>0時說明 X和 Y是正相關關係，協方差<0時 X和Y是負相關關係，協方差為0時 X和Y相互獨立。

　　協方差的計算是針對兩維的，對於n維的數據集，可以計算C(n,2)種協方差。n維數據的協方差矩陣的定義如下：　　

Dim(x)表示第x維。

對於三維(x,y,z)，其協方差矩陣如下，可看出協方差矩陣是一個對稱矩陣（symmetrical），其對角線元素為每一維的方差： 　

2.特徵向量和特徵值　

　　若

，則稱

是A的特徵值，X是對應的特徵向量。可以這樣理解：矩陣A作用在它的特徵向量X上，僅僅使得X的長度發生了變化，縮放比例就是相應的特徵值

。特徵向量只能在方陣中找到，而且並不是所有的方陣都有特徵向量，並且如果一個n*n的方陣有特徵向量，那麼就有n個特徵向量。一個矩陣的所有特徵向量是正交的，即特徵向量之間的點積為0，一般情況下，會將特徵向量歸一化，即向量長度為1。

3.PCA過程

　　第一步，獲取數據，下圖中的Data為原始數據，一共有兩個維度，可看出二維平面上的點。

　　下圖是Data在二維坐標平面上的散點圖：

　　第二步，減去平均值，對於Data中的每一維數據分別求平均值，並減去平均值，得到DataAdjust數據。

　　第三步，計算DataAdjust的協方差矩陣

　　第四步，計算協方差矩陣的特徵向量和特徵值，選取特徵向量

　　特徵值0.490833989對應的特徵向量是（-0.735178656, 0.677873399），這裡的特徵向量是正交的、歸一化的，即長度為1。

　　下圖展示DataAdjust數據和特徵向量的關係：

　　正號表示預處理後的樣本點，斜著的兩條線就分別是正交的特徵向量（由於協方差矩陣是對稱的，因此其特徵向量正交），特徵值較大的那個特徵向量是這個數據集的主要成分（principle component）。通常來說，當從協方差矩陣計算出特徵向量之後，下一步就是通過特徵值，對特徵向量進行從大到小的排序，這將給出成分意義的順序。成分的特徵值越小，其包含的信息量也就越少，因此可以適當選擇。　

　　如果數據中有n維，計算出n個特徵向量和特徵值，選擇前k個特徵向量，然後最終的數據集合只有k維，取的特徵向量命名為FeatureVector。

　　這裡特徵值只有兩個，我們選擇其中最大的那個，1.28402771，對應的特徵向量是

。

　　第五步，將樣本點投影到選取的特徵向量上，得到新的數據集

　　假設樣例數為m，特徵數為n，減去均值後的樣本矩陣為DataAdjust(m*n)，協方差矩陣是n*n，選取的k個特徵向量組成的矩陣為EigenVectors(n*k)。那麼投影后的數據FinalData為

　　這裡是FinalData(10*1) = DataAdjust(10*2矩陣)×特徵向量

　　得到結果為

　　下圖是FinalData根據最大特徵值對應的特徵向量轉化回去後的數據集形式，可看出是將DataAdjust樣本點分別往特徵向量對應的軸上做投影：

　　如果取的k=2，那麼結果是

　　可見，若使用了所有特徵向量得到的新的數據集，轉化回去之後，與原來的數據集完全一樣（只是坐標軸旋轉）。

Python實現PCA

　　將數據轉化成前K個主成分的偽碼大致如下：　

"""減去平均數計算協方差矩陣計算協方差矩陣的特徵值和特徵向量將特徵值從大到小排序保留最大的K個特徵向量將數據轉換到上述K各特徵向量構建的新空間中"""

　　代碼實現如下：

from numpy import *def loadDataSet(fileName, delim=" "): fr = open(fileName) stringArr = [line.strip().split(delim) for line in fr.readlines()] datArr = [map(float,line) for line in stringArr] return mat(datArr)def pca(dataMat, topNfeat=999999): meanVals = mean(dataMat, axis=0) DataAdjust = dataMat - meanVals #減去平均值 covMat = cov(DataAdjust, rowvar=0) eigVals,eigVects = linalg.eig(mat(covMat)) #計算特徵值和特徵向量 #print eigVals eigValInd = argsort(eigVals) eigValInd = eigValInd[:-(topNfeat+1):-1] #保留最大的前K個特徵值 redEigVects = eigVects[:,eigValInd] #對應的特徵向量 lowDDataMat = DataAdjust * redEigVects #將數據轉換到低維新空間 reconMat = (lowDDataMat * redEigVects.T) + meanVals #重構數據，用於調試 return lowDDataMat, reconMat

　　測試數據testSet.txt由1000個數據點組成。下面對數據進行降維，並用matplotlib模塊將降維後的數據和原始數據一起繪製出來。

import matplotlibimport matplotlib.pyplot as pltdataMat = loadDataSet("testSet.txt")lowDMat, reconMat = pca(dataMat,1)print "shape(lowDMat): ",shape(lowDMat)fig = plt.figure()ax = fig.add_subplot(111)ax.scatter(dataMat[:,0].flatten().A[0],dataMat[:,1].flatten().A[0],marker="^",s=90)ax.scatter(reconMat[:,0].flatten().A[0],reconMat[:,1].flatten().A[0],marker="o",s=50,c="red")plt.show()

　　結果如下圖：

Python環境

1.編譯環境：win8.1 + 32位 + python2.7

2.相關模塊安裝：

　　（1）Numpy和Scipy：numpy用於存儲和處理大型矩陣，進行數值計算。scipy是基於numpy的科學和工程計算工具，用於處理多維數組向量、矩陣、圖形（圖形圖像是像素的二維數組）、表格。下載地址：http://www.scipy.org/scipylib/download.html

　　（2）Matplotlib：python的一個圖形框架，用於數據繪圖。win32的安裝文件下載地址：http://matplotlib.org/downloads.html

　　（3）Dateutil 和 Pyparsing模塊：安裝配置matplotlib包時需要。win32的安裝文件下載地址：http://www.lfd.uci.edu/~gohlke/pythonlibs/

3.編譯遇到問題：

　　（1）提示「No module name six」，將Python27Libsite-packagesscipylib中的six.py six.pyc six.pyo三個文件拷貝到Python27Libsite-packages目錄下即可。

　　（2）提示 「ImportError:six1.3orlaterisrequired;youhave1.2.0」，說明six.py版本過低，下載最新的版本，將其中的six.py替換Python27Libsite-packages中的six.py即可，下載地址：https://pypi.python.org/pypi/six/

　　註：six模塊是Python 2和3的兼容工具

PCA可以從數據中識別其主要特徵，它是通過沿著數據最大方差方向旋轉坐標軸來實現的。其中的矩陣運算，求特徵值、特徵向量等不是很熟悉，仍需進一步學習。