利用PCA来简化数据

题外话：之前学习的时候总是喜欢手写笔记，因为喜欢自己的字（自恋脸，溜了溜了……）但是因为只有自己看所以记得乱乱的，昨天写完了人生中的第一篇博客，早上起来无意中发现居然有好几十个人浏览过，心里感觉很开心 ，以后会经常来csdn哒！
一.降维技术
先来了解下什么是降维技术，搞过机器学习的小伙伴都知道这个词，我说说自己的理解，比如说找男朋友，身高、体重、长相、学历、人品、性格、家境、户籍、职业、星座、血型等等都是姑娘们的考察方面，相信每一个姑娘对于这些方面都有自己的标准，但是如果每一方面都符合要求才行，那这男朋友太难找了，所以姑娘们只能选出自己最在意的两至三个方面重点考察，这样找到男朋友的可能性就大大提高了（这里排除即使你一无所有，我也爱你到天荒地老的那种情况），这就是所谓的降维，把之前多个维度降低到两到三个维度，简化数据，方便计算。
降维技术主要应用在数据预处理阶段，它具有以下的优点：
（a）使得数据集更易使用
（b）去除噪声，降低算法开销
（c）减少存储空间
常见的降维方法有很多种，比如主成分分析（PCA）、因子分析、独立成分分析等。今天主要学习的是主成分分析（PCA）。
二、主成分分析的算法
1.主成分分析需要解决的事
还是上面姑娘找男朋友的例子，姑娘需要选出最重要的维度，那么怎么衡量这个“最重要”呢？选择几个维度能够尽量的减少误差呢？这就是PCA需要解决的事！
2.一个简单的例子
我们先不考虑维数很大的那些高难度数据，我们先考虑如果现在有一堆二维数据，怎么把它们合理的降到一维呢？数据如下图：

如果把数据映射到u1和u2两个方向，你会选择哪一个呢？相信直观的感觉是选择u1，为什么呢？因为u1看起来更长，覆盖的长度更大，学术一点说就是样本点投射到u1上的投影的方差更大，方差其实就是在描述数据的波动程度，也可以说是数据的密集程度，那为啥方差越大越好呢？道理很简单啊，如果家里有10口人，投影到一间房子里，是20平方的房子好还是200平的房子好呢？开个玩笑啦！其实是因为区分度更大，富含的信息量更大。当然只是在此情此景下方差越大越好，其它环境下不一定。
3.[数学推导](https://www.cnblogs.com/pinard/p/6239403.html)

我的数学水平不够只能看看别人的文章，嘻嘻。大致分为这几步：
（1）数据中心化
数据中心化就是每一个数据减掉整体数据的平均值，开始我不理解为什么要这样，所以举了个例子，假如现在有4，8，9，3这四个数字，平均值是6，分别减掉平均值后变成-2，2，3，-3，好神奇！中心化之后的数字总和是0，平均值是0。
（2）标准正交基
基：在某个向量空间中，任一向量都能由一个向量组线性表示，则称这个向量组是向量空间的基。
标准正交基：这个基两两正交，且都是单位向量。
（3）新坐标系的投影方差
回想小时候学的求方差的公式，因为之前去中心化的时候得知平均值是0，所以很容易求出方差。
（4）迹
要使得投影方差的和最大，也就是最大化这个n阶矩阵的迹，这个我想了很久。首先查了什么是迹，原来迹也就是特征值的和，那么求方差的和=特征值的和，而事实上特征值确实相当于矩阵特征向量方向上的方差，所以这么算是合理的，有一遍文章特别棒，主要讲解了特征值和特征向量的意义
（5）拉格朗日函数
这个没什么，求带约束条件的函数的最值嘛，之前看SVM的时候看过更多的，更复杂的。
（6）求导
求最值问题嘛。
4.算法流程（伪代码）
前提：矩阵A是一个方阵

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

5.代码实现

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

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

def show_data():
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(dataMat[:,0].flatten().A[0],dataMat[:,1].flatten().A[0],
               marker='^', s=90)
    ax.scatter(reconMat[:,0].flatten().A[0],reconMat[:,1].flatten().A[0],
               marker='o',s=50,c='red')
    plt.show()

def pca(dataMat, topNfeat=9999999):
    meanVals = np.mean(dataMat, axis=0) #求平均值
    meanRemoved = dataMat - meanVals   #中心化
    covMat = np.cov(meanRemoved, rowvar=0) #协方差矩阵
    eigVals, eigVects = np.linalg.eig(np.mat(covMat)) #特征值和特征向量
    eigValInd = np.argsort(eigVals) #对特征值进行排序，默认从小到大
    eigValInd = eigValInd[: -(topNfeat + 1): -1] #逆序取得特征值最大的元素的索引
    redEigVects = eigVects[:, eigValInd] #用特征向量构成矩阵
    lowDDateMat = meanRemoved * redEigVects #用归一化后的各个数据与特征矩阵相乘，映射到新的空间
    reconMat = (lowDDateMat * redEigVects.T) + meanVals #还原原始数据
    return lowDDateMat, reconMat

if __name__ == "__main__":
    file = 'testSet.txt'
    dataMat = loadDataSet(file)
    lowDMat, reconMat = pca(dataMat, 1)
    print(np.shape(lowDMat))
    show_data()

输出：

(1000, 1)

可以看到降到一维的数据的直观图
6.代码分析

• split(‘\t’):按照制表符切割字符串
• strip():删除空白符（包括’\n’, ‘\r’, ‘\t’, ’ ‘)
• map(float,line):每行映射为浮点数，python3加上list，list(map(float,line))
• np.mat():将列表转换成矩阵
• flatten().A[0]
  例子：

>>> a = [[1,3],[2,4],[3,5]]
>>> a = mat(a)
>>> y = a.flatten()
>>> y
matrix([[1, 3, 2, 4, 3, 5]])
>>> y = a.flatten().A
>>> y
array([[1, 3, 2, 4, 3, 5]])
>>> shape(y)
(1, 6)
>>> shape(y[0])
(6,)
>>> y = a.flatten().A[0]
>>> y
array([1, 3, 2, 4, 3, 5])

• np.linalg.eig():求矩阵的特征值和特征向量（该函数将返回一个元组，按列排放着特征值和对应的特征向量，其中第一列为特征值，第二列为特征向量）
• np.argsort():排序函数。

 numpy.argsort(a, axis=-1, kind=’quicksort’, order=None) 
功能: 将矩阵a按照axis排序，并返回排序后的下标 
参数: a:输入矩阵， axis:需要排序的维度 
返回值: 输出排序后的下标

例子：

>>> x = np.array([[1, 5, 7], [3, 2, 4]])
>>> np.argsort(x, axis=0)
array([[0, 1, 1],
       [1, 0, 0]])  #沿着行向下(每列)的元素进行排序
>>> np.argsort(x, axis=1)
array([[0, 1, 2],
       [1, 0, 2]])  #沿着列向右(每行)的元素进行排序

• eigValInd[: -(topNfeat + 1): -1]:一开始觉得好奇，为什么非得加一个‘-’号，后来研究了一下
  例子：

   a=[[1],
       [3],
       [5],
       [4]]
    print(a[:-2:-1])
    print(a[:-4:-1])

输出：

[[4]]
[[4], [5], [3]]

7.实例（利用PCA对半导体数据降维）

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

def loadDataSet(fileName, delim='\t'):
    fr = open(fileName)
    stringArr = [line.strip().split(delim) for line in fr.readlines()]
    datArr = [list(map(float,line)) for line in stringArr]
    return np.mat(datArr)
def replaceNanWithMean():
    datMat = loadDataSet('secom.data', ' ')
    numFeat = np.shape(datMat)[1]
    for i in range(numFeat):
        meanVal = np.mean(datMat[np.nonzero(~np.isnan(datMat[:,i].A))[0],i]) #计算所有非NaN的均值
        datMat[np.nonzero(np.isnan(datMat[:,i].A))[0],i] = meanVal  #将所有NaN置为平均值
    return datMat
if __name__=="__main__":
    dataMat=replaceNanWithMean()
    meanVals = np.mean(dataMat, axis=0)
    meanRemoved = dataMat - meanVals  # 去除均值
    covMat = np.cov(meanRemoved, rowvar=0)    #计算协方差矩阵
    eigVals, eigVects = np.linalg.eig(np.mat(covMat))    #计算特征值和特征向量
    eigValInd = np.argsort(eigVals)  # 对特征值从小到大排序
    eigValInd = eigValInd[::-1]  # 反转
    sortedEigVals = eigVals[eigValInd]
    total = sum(sortedEigVals)
    varPercentage = sortedEigVals / total * 100
    eigVals_rate = 0
    eigVals_rate_list = []
    eigVals_sum = sum(sortedEigVals)
    for i in range(20):
        eigVals_rate = eigVals_rate + sortedEigVals[i] / eigVals_sum
        eigVals_rate_list.append(eigVals_rate)
    print(eigVals_rate_list)
    plt.plot(range(1, 21), eigVals_rate_list, marker='o', c='red')
    plt.show()

输出：

总结：

1.缺失值的处理

• 在590个特征下，几乎所有样本都有NaN，因此去除不完整的样本不太现实。
• 将所有的NaN替换成0，但是由于并不知道这些值的意义，所以这样做是个下策。如果它们是开氏温度，那么将它们置成0这种处理策略就太差劲了。
• 用平均值来代替缺失值。
  2.重要函数
• np.isnan(element):判断某元素是否Nan,是则返回True，不是返回False
• np.nonzero():

nonzero(a) 

返回数组a中非零元素的索引值数组。

（1）只有a中非零元素才会有索引值，那些零值元素没有索引值； 
（2）返回的索引值数组是一个2维tuple数组，该tuple数组中包含一维的array数组。其中，一维array向量的个数与a的维数是一致的。 
（3）索引值数组的每一个array均是从一个维度上来描述其索引值。比如，如果a是一个二维数组，则索引值数组有两个array，第一个array从行维度来描述索引值；第二个array从列维度来描述索引值。 
（4） 该np.transpose(np.nonzero(x)) 函数能够描述出每一个非零元素在不同维度的索引值。 
（5）通过a[nonzero(a)]得到所有a中的非零值

3.从图片可以看出前6个特征值就可以达到96%左右的方差，因此选择前6个特征值就能在一定程度上达到降维的目的。