机器学习实战（十三）——利用PCA来简化数据

机器学习实战（十三）——利用PCA来简化数据

一、降维技术

有时候数据的维度过大，导致分析时很难直观地看出关键信息；并且特征过多，很难精确提取出关键的重要特征。因此降维就显得尤为重要，将数据降低到较少的维度，有利于我们更直观地从中提取出想要的信息，同时也便于之后其他机器学习的算法进行处理。

降维技术可以应用到已标注数据和未标注数据。主要有以下几种方法：

1、主成分分析（PCA）：主要思想是将数据从原来的坐标轴转换到新的坐标轴，选择的依据主要是方差，即不断从原始数据中选取方差最大的方向，同时新坐标轴的选择要和已选择的坐标轴正交。重复次数为特征的数目。之后再选择最前面的几个新坐标轴，因为大部分的方差都在这些里面，忽略其他坐标轴，从而达到降维的效果。

在统计学习方法一书中是这样解释的：

先对数据进行标准化，之后进行正交变换，将原来由线性相关变量表示的数据通过正交变换变成若干个线性无关的新变量表示的数据。

例如上图可以这样理解：原本的数据由 $x_1{x}_2$来表示，从图中可以看出这两个特征具有一定的线性关系，或者说正相关的关系，即若已知其中一个坐标，那么对另一个坐标的预测并不是完全随机的；经过PCA处理后，由$y_1{y}_2$来表达，可以看出由这两个坐标轴表示时，数据并没有明显的相关性，或者可以说只是随机地分布在椭圆形之内，并且当已知其中一个坐标，那么对另一个坐标的预测是完全随机的。而前文所述的方差最大，其含义也可以理解成数据的随机性更强，而不是呈现一定的已知关系。

2、因子分析：假设在观察数据的生成中有一些观察不到的隐变量，即观察数据是由这些隐变量和某些噪声的线性组合，那么隐变量的维度可能会比观察数据少，就可以通过寻找隐变量来实现数据的降维

3、独立成分分析（ICA）：假设数据是由N个数据源生成的（N可以理解为要降维的目标维度），而数据是为多个数据源的混合观察结果，这些数据源之间在统计上是相关独立的。那么如果数据源的数目比较少，同样可实现降维。

二、PCA

2.1、移动坐标轴

移动坐标轴的例子在上面已经提到了。接下来来直观感受一下降维的强大力量。

例如在上图中，如果我们直接对原始数据采用某种算法来进行划分间隔，例如采用决策树（图中绿色）那么我们可以选择横坐标的特征，然后按照大小将它们划分开，但分类间隔就不会很大；其次是采用SVM（图中蓝色），可以拟合直线来区分，同时也有较大的分类间隔，但分类超平面的具体含义很难解释清楚。如果我们将坐标轴旋转：

在这个坐标轴下可以很明显地看出，我们就算采用决策树也能够有较大的分类间隔，同时其含义也非常直观。

2.2、在NumPy中实现PCA

伪代码大致如下：

去除平均值
计算协方差矩阵
计算协方差矩阵的特征值和特征向量
将特征值从大到小排序
保留最上面的N个特征向量
将数据转换到上述N个特征向量构建的新空间中

这里补充一下PCA的定理，这里x可以看成原始数据标准化后的矩阵，而y就是要转换的新的坐标轴。因此下面的代码才会涉及到对特征值排序。

具体代码为：

def loadDataSet(fileName, delim='\t'):
    fr = open(fileName)
    stringArr = [line.strip().split(delim) for line in fr.readlines()]
    dataArr = [list(map(float, line)) for line in stringArr]
    return mat(dataArr)

def PCA(dataMat, topNfeat=9999999):
    meanVals = mean(dataMat, axis=0)  # 对每一列求均值
    meanRemoved = dataMat - meanVals  # 减去均值
    covMat = cov(meanRemoved, rowvar = False)  # 求协方差矩阵，第二个参数为False代表每一列为一个变量
    eigVals,eigVects = linalg.eig(mat(covMat))  # 计算协方差矩阵的特征值和特征向量
    eigValInd = argsort(eigVals)  # 对特征值排序，但返回的是索引
    eigValInd = eigValInd[:-(topNfeat+1):-1]  # 取出我们感兴趣的前N个特征值
    redEigVects = eigVects[:,eigValInd]  # 根据索引取出前topNfeat个特征向量构成对数据转换的矩阵
    lowDDataMat = meanRemoved * redEigVects  # 计算在新空间下的坐标
    reconMat = (lowDDataMat * redEigVects.T) + meanVals  # 重构数据，降维后再返回待原先维度再加上均值
    return lowDDataMat,reconMat

解释一下这两行代码：

eigValInd = argsort(eigVals)  # 对特征值排序，但返回的是索引
eigValInd = eigValInd[:-(topNfeat+1):-1]  # 取出我们感兴趣的前N个特征值

首先，$argsort$函数将数组按照从小到大的顺序排序，返回值是索引，例子如下：

import numpy as np
x = np.array([1,4,3,-1,6,9])
x.argsort()
>>> array([3, 0, 1, 2, 4, 5], dtype=int64)

例如返回值中第一个元素为3，代表在原数组中x[3] 这个元素是最小的，然后就是 x[0] 这个元素，以此类推。

第二行代码就是取出感兴趣的前 $topNfeat$ 个最大的特征值所对应的索引位置，例子如下：

>>> A = [1,2,3,4,5,6,7,8,9]
>>> A[:3:-1]  
Out[3]: [9, 8, 7, 6, 5]
>>> A[:-2:-1]
Out[4]: [9]
>>> A[:-4:-1]
Out[5]: [9, 8, 7]
>>> A[:-1:-1]
Out[6]: []

索引的第三个数组为-1代表从后往前，因为$argsort$是从小到大；那么前面的 [: 3] 就代表从最后一个元素直到索引为3的元素，而且索引是左闭右开，索引为3这个位置是取不到的。

那么 $-topNfeat-1$ 就是从最后一个元素，取直到索引为 $-topNfeat-1$，那么就取出了最大的 $topNfeat$个特征值的索引位置。

执行完代码后，可以画个图直观展现一下降维后的数据和原始数据：

if __name__ == "__main__":
    dataMat = loadDataSet("testSet.txt")
    lowDMat,reconMat = PCA(dataMat,1)  # 两维度降低到一维度
    print(shape(lowDMat))
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(dataMat[:,0].flatten().A[0],dataMat[:,1].flatten().A[0],
               marker = '^',s = 110)
    ax.scatter(reconMat[:,0].flatten().A[0],reconMat[:,1].flatten().A[0],
               marker='o',s = 20,c = 'red')
    plt.show()

如果是降低到二维度，那么应该跟原来的数据是重合的：

三、示例：利用PCA对半导体制造数据降维

由于数据集中存在部分数据为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[nonzero(~isnan(dataMat[:,i].A))[0],i]

首先是 $isnan$，就是找出第 i 列中是Nan的那些元素的位置（是为1，不是为0）再取反就找到了不是Nan的位置，然后再索引出来求均值。

再执行如下代码即可：

if __name__  == "__main__":
	dataMat = replaceNanWithMean()
	meanVals = mean(dataMat,axis = 0)
	meanRemoved = dataMat - meanVals
	covMat = cov(meanRemoved,rowvar=False)
	eigVals,eigVects = linalg.eig(mat(covMat))
    plt.plot(eigVals[:20])  # 对前20个特征值画图
    plt.show()

可以看到：

可以看到后面那些特征值已经很小了，更有很多已经都是0。这说明只有前几个特征值能够提供更多的信息。

例如：

    # 查看前6个方差占了百分之多少
    eigVals_six = sum(eigVals[:6])
    print(eigVals_six)
    print(sum(eigVals))
    # 查看前7个方差占了百分之多少
    eigVals_six = sum(eigVals[:7])
    print(eigVals_six)
    print(sum(eigVals))
(87267225.18122156+0j)
(90146058.6410682+0j)
(87558088.73663929+0j)
(90146058.6410682+0j)

再画出占比：

	eigValsRate = []
    for i in range(30):
        temRate = float(sum(eigVals[0:i]))/ sum(eigVals)
        eigValsRate.append(temRate)
    plt.plot(eigValsRate)
    plt.show()

也可以看出前面几个特征已经使得总方差比例很接近于1了。