机器学习之PCA

1、PCA概述

在很多实际数据中，通常涉及很多的变量。大量的变量不但增加了计算的复杂度，而且有些变量有可能是噪声，
从而将数据中的主要数据“淹没”。此外虽然每一个变量都提供了相应的信息，但是很多变量可能存在一定的
相关性。因此，我们希望从数据中提取主要变量信息，用较少的新变量来表达数据中的主要信息。在主成分分析
(pricipal component analysis,PCA)中，我们使用旧变量的线性组合构建新变量。换言之，主成分分析的基本思想是利用线性变换将高维空间中的数据投影到低维空间中，并保留最主要信息。

2、PCA算法介绍

简单来说，就是将数据从原始的空间转换到新的特征空间中，例如原始的空间是三维的(x,y,z)，x,y,z分别是原始的三个基，我们可以通过某种方法，用新的坐标系(a,b,c)来表示原始的数据，那么a,b,c就是新的基，它们可以组成新的特征向量空间。在新的特征向量空间中，可能所有的数据在C上的投影都接近于0，即可以忽略，那么我们就直接用(a,b)来表示数据，这样数据就从三维的(x,y,z)降到了二维的(a,b)。

如何求新的基(a,b)?

般步骤是这样的：先对原始数据零均值化，然后求协方差矩阵，接着对协方差矩阵求特征向量和特征值，这些特征向量组成了新的特征空间。

参考：http://blog.codinglabs.org/articles/pca-tutorial.html
通俗易懂

3、选择主成分个数

我们该如何选择k，即保留多少个PCA主成分？对于高维数据来说就没那么简单：如果k过大，数据压缩率不高，在极限情况k=n 时，等于是在使用原始数据；相反地，如果k过小，那数据的近似误差太大。决定k 值时，我们通常会考虑不同 k 值可保留的方差百分比。具体来说，如果k=n ，那么我们得到的是对数据的完美近似，也就是保留了100%的方差，即原始数据的所有变化都被保留下来；相反，如果k=0，那等于是使用零向量来逼近输入数据，也就是只有0%的方差被保留下来。这个百分比也叫贡献率，通常取高于80%。

#参数：特征值，百分比
def percentage2k(featValue,percentage):
    sortArray = np.sort(featValue)    # 降序
    sortArray=sortArray[-1::-1]  #逆转，即降序 
    print("sortArray: ",sortArray)
    arraySum = sum(sortArray)
    tmpSum = 0
    num = 0
    for i in sortArray:
        tmpSum += i
        num += 1
        if tmpSum >= arraySum*percentage:
            return num

4、PCA的python实现

"""
date       2017-08
Author     ggwcr
theory     http://blog.codinglabs.org/articles/pca-tutorial.html    
Note       PCA source code 
"""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

#加载数据
def loaddata(datafile):
    return np.array(pd.read_csv(datafile,sep="\t",header=-1)).astype(np.float)


# 将数据归一化
# 返回：归一化后的矩阵与极差
def autoNorm(dataSet):
    #获取列的最小值
    minVals = dataSet.min(0)
    #获取列的最大值
    maxVals = dataSet.max(0)  
    #获得极差
    ranges = maxVals -minVals  
    #创建以传入数据集同结构的零矩阵
    normDataSet = np.zeros(np.shape(dataSet))  
    #获取列的长度
    m,n = np.shape(dataSet)
    #得到一个原矩阵与最小值相减所得的矩阵{oldValue- min )
    normDataSet = dataSet - np.tile(minVals,(m,1))
    #计算{oldValue- min ) / (max-min)
    normDataSet = normDataSet/np.tile(ranges,(m,1))
    return normDataSet,ranges

# 零均值化
def zeroMean(dataSet):
    meanVal = np.mean(dataSet,axis=0)  # axis = 0 表示按照列来求均值 
    ZMDataSet = dataSet - meanVal
    return ZMDataSet,meanVal

# 指定k值
def pca(dataSet,k):
    # 先将数据归一化
    normDataSet,ranges = autoNorm(dataSet)
    # 对归一化后的数据进行零均值化
    ZMDataSet,meanVal = zeroMean(normDataSet)
    m,n = np.shape(ZMDataSet)
    # 求零均值化后的协方差矩阵
    covMat = np.cov(ZMDataSet,rowvar=0)  #若rowvar=0，说明传入的数据一行代表一个样本
    #求解协方差矩阵的特征值和特征向量
    #返回值为：特征值和特征向量
    featValue, featVec=  np.linalg.eig(np.mat(covMat))  
    #按照featValue进行从大到小排序
    index = np.argsort(-featValue) 
    finalData = []
    # 对需要降到的维数大小进行判断
    if k > n:
        print ("kmust lower than feature number")
        return
    else:
        #注意特征向量时列向量，而numpy的二维矩阵(数组)a[m][n]中，a[1]表示第1行值
        selectVec = np.matrix(featVec.T[index[:k]]) #所以这里需要进行转置
        lowDataMat = ZMDataSet* selectVec.T 
        reconDataMat = (lowDataMat* selectVec) +meanVal  
    return lowDataMat, reconDataMat    

# 根据贡献率，自己找k值
def pcaNK(dataSet,percentage):
    # 先将数据归一化
    normDataSet,ranges = autoNorm(dataSet)
    # 对归一化后的数据进行零均值化
    ZMDataSet,meanVal = zeroMean(normDataSet)
    m,n = np.shape(ZMDataSet)
    # 求零均值化后的协方差矩阵
    covMat = np.cov(ZMDataSet,rowvar=0)  #若rowvar=0，说明传入的数据一行代表一个样本
    #求解协方差矩阵的特征值和特征向量
    #返回值为：特征值和特征向量
    featValue, featVec=  np.linalg.eig(np.mat(covMat))  
    print("featValue: ",featValue)
    #根据percentage求k值
    k = percentage2k(featValue,percentage)
    print("k is :",k)
    #按照featValue进行从大到小排序
    index = np.argsort(-featValue) 
    finalData = []
    #注意特征向量时列向量，而numpy的二维矩阵(数组)a[m][n]中，a[1]表示第1行值
    selectVec = np.matrix(featVec.T[index[:k]]) #所以这里需要进行转置
    lowDataMat = ZMDataSet* selectVec.T 
    reconDataMat = (lowDataMat* selectVec) +meanVal  
    return lowDataMat, reconDataMat 

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() 

def test():
    X = [[2.5, 0.5, 2.2, 1.9, 3.1, 2.3, 2, 1, 1.5, 1.1],
        [2.4, 0.7, 2.9, 2.2, 3.0, 2.7, 1.6, 1.1, 1.6, 0.9]]
    XMat = np.matrix(X).T  
    k = 2
    return pca(XMat, k)

#根据数据集data.txt
def main():    
    datafile = "c:\data.csv"
    XMat = loaddata(datafile)
    k = 2
    return pca(XMat, k)
if __name__ == "__main__":
    finalData, reconMat = main()
    plotBestFit(finalData, reconMat)

蓝色部分为重构后的原始数据，红色则是提取后的二维特征！