PCA主成分分析Python實現

作者：拾毅者

出處： http://www.csuldw.com/2016/01/22/2016-02-28-pca/

</div>

PCA（principle component analysis） ，主成分分析，主要是用來降低數據集的維度，然后挑選出主要的特征。原理簡單，實現也簡單。關于原理公式的推導，本文不會涉及，因為有很多的大牛都已經謝過了，這里主要關注下實現，算是鍛煉一下自己。

本來是在復習LDA的，然后就看到了PCA，就跟著下面這篇文章的步驟，把PCA用python實現了一遍，具體的思想可以參考這篇文章，講的通俗易懂，主要是有個實例參考，值得擁有！

• JerryLead之PCA主成分分析

下面自己來簡單的清理下思路！

PCA思想

移動坐標軸，將n維特征映射到k維上（k<n），這k維是全新的正交特征。這k維特征稱為主元，是重新構造出來的k維特征，而不是簡單地從n維特征中去除其余n-k維特征。

說到PCA難免會提到LDA（linear discriminate analysis，線性判別分析），以及FA（factor analysis，因子分析）。關于LDA，打算有時間也用代碼實現一遍，下面給出它的主要思想。

LDA思想：最大類間距離，最小類內距離。簡而言之，第一，為了實現投影后的兩個類別的距離較遠，用映射后兩個類別的均值差的絕對值來度量。第二，為了實現投影后，每個類內部數據點比較聚集，用投影后每個類別的方差來度量。

三者的描述如下

以下內容引自 Wikipedia- Linear discriminant analysis

區別：PCA選擇樣本點投影具有最大方差的方向，LDA選擇分類性能最好的方向。

好了，下面來看下實現源碼！

基本步驟

基本步驟：

• 對數據進行歸一化處理（代碼中并非這么做的，而是直接減去均值）
• 計算歸一化后的數據集的協方差矩陣
• 計算協方差矩陣的特征值和特征向量
• 保留最重要的k個特征（通常k<n），可以自己制定，也可以選擇個閾值，讓后通過前k個特征值之和減去后面n-k個特征值之和大于這個閾值，找到這個k
• 找出k個特征值對應的特征向量
• 將m n的數據集乘以k個n維的特征向量的特征向量（n k）,得到最后降維的數據。

源碼實現

1.首先引入numpy，由于測試中用到了pandas和matplotlib，所以這里一并加載

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

#計算均值,要求輸入數據為numpy的矩陣格式，行表示樣本數，列表示特征 
def meanX(dataX):
    return np.mean(dataX,axis=0)#axis=0表示按照列來求均值，如果輸入list,則axis=1

"""
參數：
 - XMat：傳入的是一個numpy的矩陣格式，行表示樣本數，列表示特征 
 - k：表示取前k個特征值對應的特征向量
返回值：
 - finalData：參數一指的是返回的低維矩陣，對應于輸入參數二
 - reconData：參數二對應的是移動坐標軸后的矩陣
"""
def pca(XMat, k):
    average = meanX(XMat) 
    m, n = np.shape(XMat)
    data_adjust = []
    avgs = np.tile(average, (m, 1))
    data_adjust = XMat - avgs
    covX = np.cov(data_adjust.T)   #計算協方差矩陣
    featValue, featVec=  np.linalg.eig(covX)  #求解協方差矩陣的特征值和特征向量
    index = np.argsort(-featValue) #按照featValue進行從大到小排序
    finalData = []
    if k > n:
        print "k must lower than feature number"
        return
    else:
        #注意特征向量時列向量，而numpy的二維矩陣(數組)a[m][n]中，a[1]表示第1行值
        selectVec = np.matrix(featVec.T[index[:k]]) #所以這里需要進行轉置
        finalData = data_adjust * selectVec.T 
        reconData = (finalData * selectVec) + average  
    return finalData, reconData

#輸入文件的每行數據都以\t隔開
def loaddata(datafile):
    return np.array(pd.read_csv(datafile,sep="\t",header=-1)).astype(np.float)

def plotBestFit(data1, data2):    
    dataArr1 = np.array(data1)
    dataArr2 = np.array(data2)

    m = np.shape(dataArr1)[0]
    axis_x1 = []
    axis_y1 = []
    axis_x2 = []
    axis_y2 = []
    for i in range(m):
        axis_x1.append(dataArr1[i,0])
        axis_y1.append(dataArr1[i,1])
        axis_x2.append(dataArr2[i,0]) 
        axis_y2.append(dataArr2[i,1])                 
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(axis_x1, axis_y1, s=50, c='red', marker='s')
    ax.scatter(axis_x2, axis_y2, s=50, c='blue')
    plt.xlabel('x1'); plt.ylabel('x2');
    plt.savefig("outfile.png")
    plt.show()

#根據數據集data.txt
def main():    
    datafile = "data.txt"
    XMat = loaddata(datafile)
    k = 2
    return pca(XMat, k)
if __name__ == "__main__":
    finalData, reconMat = main()
    plotBestFit(finalData, reconMat)