Python 主成分分析PCA

　　主成分分析（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='\t'):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 eigValseigValInd = 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 matplotlib
import 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”，將\Python27\Lib\site-packages\scipy\lib中的six.py six.pyc six.pyo三個文件拷貝到\Python27\Lib\site-packages目錄下即可。

　　（2）提示 “ImportError:?six?1.3?or?later?is?required;?you?have?1.2.0”，說明six.py版本過低，下載最新的版本，將其中的six.py替換\Python27\Lib\site-packages中的six.py即可，下載地址：https://pypi.python.org/pypi/six/

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

?

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