機器學習實戰——PCA（主成分分析）
本章關於PCA的程式碼雖少，但涉及到的知識卻很多，由於數學知識比較淺薄，所以在看這章時提前查詢資料複習了很多的概率論和統計學知識和python基礎知識，這裡記錄的很多都是關於PCA的相關知識或理論（例如：特徵向量、協方差矩陣等），由於部分知識涉及較多，講的有點詳細所以文章篇幅較長儘量縮減了，下面進入正文。
通常我們可以很清楚的看到一維資料，或直觀的二維圖形，但實際中很多的資料遠不止1、2、3維，這些資料往往擁有超出顯示能力的更多特徵。資料顯示並非大規模特徵下的唯一難題，對資料進行簡化還有如下一系列原因：
- 使得資料集更易使用；
- 降低很多演算法的計算開銷；
- 去除噪聲點；
- 使得結果易懂。
降維技術：
主成分分析（Principal Component Analysis, PCA）。資料從原來的座標系轉化到了新的座標系，新座標系的選擇是由資料本身決定的。第一個新座標軸選擇的是原始資料中方差最大的方向，第二個新座標軸的選擇和第一個座標軸正交且具有最大方差的方向。重複選擇新座標軸，最終大部分方差都包含在最前面的幾個新座標軸中。。因此，忽略餘下的座標軸，達到降維的目的。
因子分析（Factor Analysis）。在因子分析中，我們假定在觀察資料的生成中有一些觀察不到的隱變數（latent variable）。假設觀察資料時這些隱變數和某些噪聲的線性組合。通過找到隱變數就可以實現資料的降維。
獨立主成分分析（Independent Component Analysis, ICA）。ICA假設資料是從N個數據源生成的。假設資料位多個數據源的混合觀察結果，這些資料來源之間在統計上是相互獨立的，而在PCA中假設資料是不相關的。同因子分析一樣，如果資料來源的數目少於觀察資料的數目，則可以實現降維過程。
流行學習（manifold learning）。
在這些流行的降維技術中，PCA是使用的最廣泛的資料壓縮演算法，它的原理上面已經提到，就是移動座標軸，而移動是由資料本身決定的（即每次都是尋找資料最大方差方向），降維後我們可以使用更簡單的分類器或者得到更好的分類間隔，下圖說明：
這些樣本的協方差矩陣
此時的三個特徵向量相互正交，構成了三維空間下的座標（v1、v2、v3）。
討論了通過座標軸旋轉的方法，但是這並沒有解決降維的問題，如何通過座標軸的旋轉來降維呢？
前面提到的第一個主成分就是從資料差異性最大的（即方差最大）的方向提取出來的，第二個主成分則來自於資料差異性次大的方向，並且與第一個主成分方向正交。通過資料集的協方差矩陣及其特徵值分析，我們就可以求得這些主成分的值。
一旦得到了協方差矩陣的特徵向量，我們就可以保留最大的N個值。這些特徵向量也給出了N個最重要特徵的真實結構。我麼可以通過將資料乘上這N個特徵向量而將它轉換到新的空間，以此達到了降維的目的。
補充相關知識：
好了上面已經把PCA降維的原理說了一遍，下面補充相關知識更仔細的解釋請參考這裡（如果知道了可以跳過）：
方差：
最先提到的一個概念，也是旋轉座標軸的依據。之所以使用方差作為旋轉條件是因為：最大方差給出了資料的最重要的資訊。
協方差：
用來衡量兩個變數的總體誤差，方差是協方差的一種特殊情況，即當兩個變數相同。可以通俗的理解為：兩個變數在變化過程中是否同向變化？還是反方向變化？同向或反向程度如何？取值為負∞到正∞仿照方差的定義，度量各個維度偏離其均值的程度，定義為：
由協方差的定義可以推出兩個性質：
協方差矩陣：
協方差只能處理二維問題（即兩個特徵X、Y），維數多了自然需要計算多個協方差，比如n維的資料集需要計算 個協方差，自然而然我們會想到使用矩陣來組織這些資料，協方差定義：
可見，協方差是一個對稱的矩陣，且對角線是各維度上的方差。正是由於協方差矩陣為對稱矩陣所以矩陣分解後特徵值所對應的特徵向量一定無線性關係，且相互之間一定正交，即內積為零。
特徵向量：
對協方差矩陣的特徵向量最直觀的解釋之一是它總是指向資料方差最大的方向（上面的u、v）。所以我們需要求得協方差矩陣，然後計算出其特徵向量，通過對特徵值的排序，選出我們要求的N個特徵向量（即N個最重要特徵的真實結構），用原資料乘上這N個特徵向量而將它轉換到新的空間中。（在numpy中linalg的eig方法可以求得特徵值、特徵向量，對特徵值排序後選擇最大的特徵向量）
將二維特徵降轉為一維效果如下圖（其它維數腦補）：
python程式碼實現PCA降維處理：
虛擬碼：
- 去除平均值
- 計算協方差矩陣
- 計算協方差矩陣的特徵值和特徵向量
- 將特徵值從大到小排序
- 保留最上面的N個特徵向量
- 將資料轉換到上述N個特徵向量構建的新空間中
程式碼：
樣例資料：
from numpy import *
import matplotlib.pyplot as plt
# 載入本地資料
def loadDataSet(fileName,delim='\t'):
fr = open(fileName)
stringArr = [line.strip().split(delim) for line in fr.readlines()]
datArr = [list(map(float,line)) for line in stringArr] # py3需要注意使用list改變資料型別
return mat(datArr) # 轉換為二維的陣列便於運算
# pca降維函式 （原始資料 壓縮維數） 不指定壓縮維數預設9999999即維數全部返回
def pca(dataMat,topNfeat=9999999):
meanVals = mean(dataMat,axis=0) # 按列求均值
meanRemoved = dataMat - meanVals # 去平均值每個值減去該列的平均值（後面有解釋為什麼要去均值）
#計算協方差矩陣，除數n-1是為了得到協方差的無偏估計
#cov(X,0) = cov(X) 除數是n-1(n為樣本個數)
#cov(X,1) 除數是n
covMat = cov(meanRemoved,rowvar=0) # 協方差矩陣計算 （rowvar=0 後面有解釋這個引數的重要意義）
eigVals,eigVects = linalg.eig(mat(covMat)) # 分解計算協方差矩陣的（特徵值 特徵向量）
eigValInd = argsort(eigVals) # 獲取從小到大排序特徵值的下標 於sort有差別
eigValInd = eigValInd[:-(topNfeat+1):-1] # 獲得最大的特徵值的下標N個 (起始位置:結束位置:步長 步長為負則從右往左取；結束位置步長同時為負則從右排列向左排列)
redEigVects = eigVects[:,eigValInd] # 保留N個特徵值最大列的特徵向量 (即資料方差最大的N個方向)
lowDDataMat = meanRemoved * redEigVects # 將資料轉換到新空間 （用特徵向量*資料集）
reconMat = (lowDDataMat * redEigVects.T) + meanVals # 利用降維後的矩陣反構出原始矩陣（作用：用作測試，可以和未壓縮的原矩陣比較）
return lowDDataMat,reconMat # 壓縮後的矩陣 反構出的原始矩陣
dataMat = loadDataSet('testSet.txt')
lowDMat,reconMat = pca(dataMat,1) # 返回了 降維後的資料 反構出的原始資料
# 畫圖
fig = plt.figure()
ax = fig.add_subplot(111)
# ax.scatter(lowDMat[:,0].flatten().A[0],zeros(len(lowDMat))) # 降維後的圖
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()解釋：減去均值後再除以標準差得出的數值就是標準化資料，例如：影響某國居民體重的兩個因素是收入和運動量，這兩個因素的均值分別是50000元和50小時，標準差是2000元和1小時，那麼你就很難比較這兩因素的變動對體重的影響，所以就要標準差，減去均值再除以標準差得出的資料就表示“偏離均值多少個標準差”，就能夠在迴歸時比較各自對體重的影響程度！
三角形：原始資料 圓形：第一主成分
降維資料如下圖：
示例：利用PCA對半導體制造資料降維（590維特徵）
（呼叫函式之前需要構造一個函式來處理資料中的空值NaN，這裡使用均值代替。我試過pandas資料來處理很方便也很好理解，但是執行速度慢了太多。）
# 示例 半導體資料降維
# 將NaN值置換為平均值
def replaceNanWithMean():
datMat = loadDataSet('secom.data',' ') # 解析資料
numFeat = shape(datMat)[1] # 獲取特徵維度
for i in range(numFeat):
meanVal = mean(datMat[nonzero(~isnan(datMat[:,i].A))[0],i]) # 利用該維度所有非NaN特徵求取均值
datMat[nonzero(isnan(datMat[:,i].A))[0],i] = meanVal # 將該維度中所有NaN特徵全部用均值替換
return datMat
dataMat = replaceNanWithMean()
meanVals = mean(dataMat, axis=0)
meanRemoved = dataMat - meanVals
covMat = cov(meanRemoved, rowvar=0)
eigVals, eigVects = linalg.eig(mat(covMat))
print(sum(eigVals)*0.9) # 計算90%的主成分方差總和
print(sum(eigVals[:6])) # 計算前6個主成分所佔的方差
plt.plot(eigVals[:20]) # 對前20個畫圖觀察
plt.show()
特徵值中超過20%的值都為0：
sklearn的PCA處理：
同時我們使用sklearn中PCA處理方法來做步驟更為簡單的處理，並與上面的結果作比較：
from sklearn.decomposition import PCA
pca = PCA() # n_components引數選擇降維程度
pca = pca.fit(replaceNanWithMean()) # fit_transform()為轉換資料
main_var = pca.explained_variance_ # 特徵值
print(sum(main_var)*0.9) #下面與手寫程式碼作比較
print(sum(main_var[:6]))
plt.plot(main_var[:20])
plt.show()
將兩種方法作比較：
效果不錯，與上面作者提供的程式碼相差無幾,前6個主成分覆蓋了96.8%的方差，前20個主成分覆蓋了99.3%的方差。如果保留前6個而去除後584個主成分，就可以實現大概100：1的壓縮比。有效的主成分數目取決於資料集和具體應用。上面顯示我們編寫的程式得到主成分的90%為81131452.77696137與sklearn計算得出的81131452.7769614幾乎一致。前6個也幾乎一樣。
以上為全部內容，部分知識為個人理解，如有問題，請指出，在學習過程中參考了一些資料，十分感謝，如果由需要可以前去檢視。
主要參考書籍：《機器學習實戰》