(3)主成分分析(PCA)——基於python+numpy
阿新 • • 發佈:2018-12-31
【機器學習演算法實現】主成分分析(PCA)——基於python+numpy
1、PCA演算法介紹
主成分分析(Principal Components Analysis),簡稱PCA,是一種資料降維技術,用於資料預處理。一般我們獲取的原始資料維度都很高,比如1000個特徵,在這1000個特徵中可能包含了很多無用的資訊或者噪聲,真正有用的特徵才100個,那麼我們可以運用PCA演算法將1000個特徵降到100個特徵。這樣不僅可以去除無用的噪聲,還能減少很大的計算量。
PCA演算法是如何實現的?
簡單來說,就是將資料從原始的空間中轉換到新的特徵空間中,例如原始的空間是三維的(x,y,z),x、y、z分別是原始空間的三個基,我們可以通過某種方法,用新的座標系(a,b,c)來表示原始的資料,那麼a、b、c就是新的基,它們組成新的特徵空間。在新的特徵空間中,可能所有的資料在c上的投影都接近於0,即可以忽略,那麼我們就可以直接用(a,b)來表示資料,這樣資料就從三維的(x,y,z)降到了二維的(a,b)。
問題是如何求新的基(a,b,c)?
一般步驟是這樣的:先對原始資料零均值化,然後求協方差矩陣,接著對協方差矩陣求特徵向量和特徵值,這些特徵向量組成了新的特徵空間。具體的細節,推薦Andrew Ng的網頁教程:Ufldl 主成分分析 ,寫得很詳細。
2、PCA演算法實現
語言:Python 函式庫:Numpy- >>> import numpy as np
根據上面提到的一般步驟來實現PCA演算法
(1)零均值化
假如原始資料集為矩陣dataMat,dataMat中每一行代表一個樣本,每一列代表同一個特徵。零均值化就是求每一列的平均值,然後該列上的所有數都減去這個均值。也就是說,這裡零均值化是對每一個特徵而言的,零均值化都,每個特徵的均值變成0。實現程式碼如下:- def zeroMean(dataMat):
- meanVal=np.mean(dataMat,axis=0) #按列求均值,即求各個特徵的均值
- newData=dataMat-meanVal
- return newData,meanVal
(2)求協方差矩陣
- newData,meanVal=zeroMean(dataMat)
-
covMat=np.cov(newData,rowvar=0
numpy中的cov函式用於求協方差矩陣,引數rowvar很重要!若rowvar=0,說明傳入的資料一行代表一個樣本,若非0,說明傳入的資料一列代表一個樣本。因為newData每一行代表一個樣本,所以將rowvar設定為0。 covMat即所求的協方差矩陣。
(3)求特徵值、特徵矩陣
呼叫numpy中的線性代數模組linalg中的eig函式,可以直接由covMat求得特徵值和特徵向量:- eigVals,eigVects=np.linalg.eig(np.mat(covMat))
(4)保留主要的成分[即保留值比較大的前n個特徵]
第三步得到了特徵值向量eigVals,假設裡面有m個特徵值,我們可以對其排序,排在前面的n個特徵值所對應的特徵向量就是我們要保留的,它們組成了新的特徵空間的一組基n_eigVect。將零均值化後的資料乘以n_eigVect就可以得到降維後的資料。程式碼如下:- eigValIndice=np.argsort(eigVals) #對特徵值從小到大排序
- n_eigValIndice=eigValIndice[-1:-(n+1):-1] #最大的n個特徵值的下標
- n_eigVect=eigVects[:,n_eigValIndice] #最大的n個特徵值對應的特徵向量
- lowDDataMat=newData*n_eigVect #低維特徵空間的資料
- reconMat=(lowDDataMat*n_eigVect.T)+meanVal #重構資料
- return lowDDataMat,reconMat
程式碼中有幾點要說明一下,首先argsort對特徵值是從小到大排序的,那麼最大的n個特徵值就排在後面,所以eigValIndice[-1:-(n+1):-1]就取出這個n個特徵值對應的下標。【python裡面,list[a:b:c]代表從下標a開始到b,步長為c。】 reconMat是重構的資料,乘以n_eigVect的轉置矩陣,再加上均值meanVal。 OK,這四步下來就可以從高維的資料dataMat得到低維的資料lowDDataMat,另外,程式也返回了重構資料reconMat,有些時候reconMat課便於資料分析。 貼一下總的程式碼:
- #零均值化
- def zeroMean(dataMat):
- meanVal=np.mean(dataMat,axis=0) #按列求均值,即求各個特徵的均值
- newData=dataMat-meanVal
- return newData,meanVal
- def pca(dataMat,n):
- newData,meanVal=zeroMean(dataMat)
- covMat=np.cov(newData,rowvar=0) #求協方差矩陣,return ndarray;若rowvar非0,一列代表一個樣本,為0,一行代表一個樣本
- eigVals,eigVects=np.linalg.eig(np.mat(covMat))#求特徵值和特徵向量,特徵向量是按列放的,即一列代表一個特徵向量
- eigValIndice=np.argsort(eigVals) #對特徵值從小到大排序
- n_eigValIndice=eigValIndice[-1:-(n+1):-1] #最大的n個特徵值的下標
- n_eigVect=eigVects[:,n_eigValIndice] #最大的n個特徵值對應的特徵向量
- lowDDataMat=newData*n_eigVect #低維特徵空間的資料
- reconMat=(lowDDataMat*n_eigVect.T)+meanVal #重構資料
- return lowDDataMat,reconMat
3、選擇主成分個數
文章寫到這裡還沒有完,應用PCA的時候,對於一個1000維的資料,我們怎麼知道要降到幾維的資料才是合理的?即n要取多少,才能保留最多資訊同時去除最多的噪聲?一般,我們是通過方差百分比來確定n的,這一點在Ufldl教程中說得很清楚,並且有一條簡單的公式,下面是該公式的截圖:根據這條公式,可以寫個函式,函式傳入的引數是百分比percentage和特徵值向量,然後根據percentage確定n,程式碼如下:
- def percentage2n(eigVals,percentage):
- sortArray=np.sort(eigVals) #升序
- sortArray=sortArray[-1::-1] #逆轉,即降序
- arraySum=sum(sortArray)
- tmpSum=0
- num=0
- for i in sortArray:
- tmpSum+=i
- num+=1
- if tmpSum>=arraySum*percentage:
- return num
那麼pca函式也可以重寫成百分比版本,預設百分比99%。
- def pca(dataMat,percentage=0.99):
- newData,meanVal=zeroMean(dataMat)
- covMat=np.cov(newData,rowvar=0) #求協方差矩陣,return ndarray;若rowvar非0,一列代表一個樣本,為0,一行代表一個樣本
- eigVals,eigVects=np.linalg.eig(np.mat(covMat))#求特徵值和特徵向量,特徵向量是按列放的,即一列代表一個特徵向量
- n=percentage2n(eigVals,percentage) #要達到percent的方差百分比,需要前n個特徵向量
- eigValIndice=np.argsort(eigVals) #對特徵值從小到大排序
- n_eigValIndice=eigValIndice[-1:-(n+1):-1] #最大的n個特徵值的下標
- n_eigVect=eigVects[:,n_eigValIndice] #最大的n個特徵值對應的特徵向量
- lowDDataMat=newData*n_eigVect #低維特徵空間的資料
- reconMat=(lowDDataMat*n_eigVect.T)+meanVal #重構資料
- return lowDDataMat,reconMat