三種方法實現(xiàn)PCA算法(Python)

??主成分分析,即Principal Component Analysis(PCA),是多元統(tǒng)計中的重要內(nèi)容,也廣泛應(yīng)用于機器學(xué)習(xí)和其它領(lǐng)域。它的主要作用是對高維數(shù)據(jù)進行降維。PCA把原先的n個特征用數(shù)目更少的k個特征取代,新特征是舊特征的線性組合,這些線性組合最大化樣本方差,盡量使新的k個特征互不相關(guān)。關(guān)于PCA的更多介紹,請參考:https://en.wikipedia.org/wiki/Principal_component_analysis.
??PCA的主要算法如下:

  • 組織數(shù)據(jù)形式,以便于模型使用;
  • 計算樣本每個特征的平均值;
  • 每個樣本數(shù)據(jù)減去該特征的平均值(歸一化處理);
  • 求協(xié)方差矩陣;
  • 找到協(xié)方差矩陣的特征值和特征向量;
  • 對特征值和特征向量重新排列(特征值從大到小排列);
  • 對特征值求取累計貢獻率;
  • 對累計貢獻率按照某個特定比例,選取特征向量集的字跡合;
  • 對原始數(shù)據(jù)(第三步后)。

??其中協(xié)方差矩陣的分解可以通過按對稱矩陣的特征向量來,也可以通過分解矩陣的SVD來實現(xiàn),而在Scikit-learn中,也是采用SVD來實現(xiàn)PCA算法的。關(guān)于SVD的介紹及其原理,可以參考:矩陣的奇異值分解(SVD)(理論)。
??本文將用三種方法來實現(xiàn)PCA算法,一種是原始算法,即上面所描述的算法過程,具體的計算方法和過程,可以參考:A tutorial on Principal Components Analysis, Lindsay I Smith. 一種是帶SVD的原始算法,在Python的Numpy模塊中已經(jīng)實現(xiàn)了SVD算法,并且將特征值從大從小排列,省去了對特征值和特征向量重新排列這一步。最后一種方法是用Python的Scikit-learn模塊實現(xiàn)的PCA類直接進行計算,來驗證前面兩種方法的正確性。
??用以上三種方法來實現(xiàn)PCA的完整的Python如下:

import numpy as np
from sklearn.decomposition import PCA
import sys
#returns choosing how many main factors
def index_lst(lst, component=0, rate=0):
    #component: numbers of main factors
    #rate: rate of sum(main factors)/sum(all factors)
    #rate range suggest: (0.8,1)
    #if you choose rate parameter, return index = 0 or less than len(lst)
    if component and rate:
        print('Component and rate must choose only one!')
        sys.exit(0)
    if not component and not rate:
        print('Invalid parameter for numbers of components!')
        sys.exit(0)
    elif component:
        print('Choosing by component, components are %s......'%component)
        return component
    else:
        print('Choosing by rate, rate is %s ......'%rate)
        for i in range(1, len(lst)):
            if sum(lst[:i])/sum(lst) >= rate:
                return i
        return 0

def main():
    # test data
    mat = [[-1,-1,0,2,1],[2,0,0,-1,-1],[2,0,1,1,0]]
    
    # simple transform of test data
    Mat = np.array(mat, dtype='float64')
    print('Before PCA transforMation, data is:\n', Mat)
    print('\nMethod 1: PCA by original algorithm:')
    p,n = np.shape(Mat) # shape of Mat 
    t = np.mean(Mat, 0) # mean of each column
    
    # substract the mean of each column
    for i in range(p):
        for j in range(n):
            Mat[i,j] = float(Mat[i,j]-t[j])
            
    # covariance Matrix
    cov_Mat = np.dot(Mat.T, Mat)/(p-1)
    
    # PCA by original algorithm
    # eigvalues and eigenvectors of covariance Matrix with eigvalues descending
    U,V = np.linalg.eigh(cov_Mat) 
    # Rearrange the eigenvectors and eigenvalues
    U = U[::-1]
    for i in range(n):
        V[i,:] = V[i,:][::-1]
    # choose eigenvalue by component or rate, not both of them euqal to 0
    Index = index_lst(U, component=2)  # choose how many main factors
    if Index:
        v = V[:,:Index]  # subset of Unitary matrix
    else:  # improper rate choice may return Index=0
        print('Invalid rate choice.\nPlease adjust the rate.')
        print('Rate distribute follows:')
        print([sum(U[:i])/sum(U) for i in range(1, len(U)+1)])
        sys.exit(0)
    # data transformation
    T1 = np.dot(Mat, v)
    # print the transformed data
    print('We choose %d main factors.'%Index)
    print('After PCA transformation, data becomes:\n',T1)
    
    # PCA by original algorithm using SVD
    print('\nMethod 2: PCA by original algorithm using SVD:')
    # u: Unitary matrix,  eigenvectors in columns 
    # d: list of the singular values, sorted in descending order
    u,d,v = np.linalg.svd(cov_Mat)
    Index = index_lst(d, rate=0.95)  # choose how many main factors
    T2 = np.dot(Mat, u[:,:Index])  # transformed data
    print('We choose %d main factors.'%Index)
    print('After PCA transformation, data becomes:\n',T2)
    
    # PCA by Scikit-learn
    pca = PCA(n_components=2) # n_components can be integer or float in (0,1)
    pca.fit(mat)  # fit the model
    print('\nMethod 3: PCA by Scikit-learn:')
    print('After PCA transformation, data becomes:')
    print(pca.fit_transform(mat))  # transformed data
            
main()

運行以上代碼,輸出結(jié)果為:

Eclipse運行結(jié)果

??這說明用以上三種方法來實現(xiàn)PCA都是可行的。這樣我們就能理解PCA的具體實現(xiàn)過程啦~~
??有興趣的讀者可以用其它語言實現(xiàn)一下哈。


參考文獻:

  1. PCA 維基百科: https://en.wikipedia.org/wiki/Principal_component_analysis.
  2. 講解詳細又全面的PCA教程: A tutorial on Principal Components Analysis, Lindsay I Smith.
  3. 博客:矩陣的奇異值分解(SVD)(理論):http://www.cnblogs.com/jclian91/p/8022426.html.
  4. 博客:主成分分析PCA: https://www.cnblogs.com/zhangchaoyang/articles/2222048.html.
  5. Scikit-learn的PCA介紹:http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html.
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容