基于Python和Numpy的主成分分析

参考文献：

http://www.cnblogs.com/jerrylead/archive/2011/04/18/2020209.html

https://my.oschina.NET/gujianhan/blog/225241

http://blog.csdn.net/jinshengtao/article/details/18599165

http://blog.csdn.Net/u012162613/article/details/42177327

http://deeplearning./wiki/index.PHP/%E4%B8%BB%E6%88%90%E5%88%86%E5%88%86%E6%9E%90

主成分分析（Principal Component Analysis）是一种对特征进行降维的方法。在机器学习以及影像分类中广泛使用，能够极大提高无监督特征学习的速度。由于观测指标间存在相关性，将导致信息的重叠与低效。为此，我们更倾向于用更少量的、尽可能多的能反映源特征的新特征来代替原始观测指标。这就是主成分分析，其可以看作是高维空间通过旋转坐标系，找到最佳投影的过程。换言之，即将n维特征映射到k维上（k<n），k维特征成为主成分，其包含了n维特征的绝大部分性质，去除了重复性质。

从实际算法实习上来看，PCA主要分为三个部分。（1）生成协方差矩阵；（2）计算特征值和特征向量，并选取主成分；（3）将原始数据投影到降维的子空间中。本文主要通过从对图像进行PCA处理的思路展开。

（1） 第一步生成协方差矩阵首先，什么是协方差矩阵？首先说方差，当我们衡量一组数据的离散程度时，使用方差来表示。即如下图所示。

S为方差。即样本中各个数据与其平均值之差的平方的和的平方。在matlab或者numpy中可以利用cov（X,X）计算。但是，该方法只能描述数据自身的离散程度。换句话说，只能描述一维数据。当我们想研究多维数据之间的关系时，就会用到协方差。如下图所示。

当我们研究维数大于2的数据组之间的关系时，便需要用到协方差矩阵。如C表示3维数据的协方差矩阵，对角线上为X,Y,Z各自的方法，其他位置表示数据之间的协方差。协方差越小，数据越相关。

那么如何计算协方差矩阵，matlab和numpy都可以利用COV（X）进行直接计算。注意这个地方输入的X为一个矩阵，在matlab中默认每一列为一个一维数据，行数代表了数据组的维数。其实现方法主要有以下两种：一是依据定义计算，即将每一组数据减去该数据的平均值再与另一组数据的结果相乘，除以n-1.第二种是简化的去中心化方法。值得注意的是numpy中的cov函数与matlab不同，其将每一行作为一个一维数据。因此利用cov进行计算，需先对其转置。其代码如下。

[python] view plain copy

1. #-*- coding:utf-8 -*-  
2. from numpy import *  
3. I=array([[1,2,3],[2,3,3],[2,2,2]])  
4. #依据定义计算  
5. m,n = I.shape  
6. I1 = zeros(I.shape)  
7. for i in range(0, n):  
8.     for j in range(0, n):  
9.         I1[i, j] = sum(dot((I[:, i] - mean(I[:, i])), (I[:, j] - mean(I[:, j])))) / (m - 1)  
10. #简化去中心化方法  
11. mean_n=I.mean(axis=0)#求解输入数据的平均值  
12. Im=I-mean_n#数据中心化  
13. M=(dot(Im.T,Im))/(m-1)#计算协方差矩阵  
14. #cov函数  
15. I2=cov(I.T)  
16. print I1,'\n',M,'\n',I2  

（2） 第二步是计算特征值和特征向量，并选取主成分。

Numpy中计算特征值和特征向量的函数主要有两个，eig（）和eigh（）。一个计算的特征值为复数形式，一个为实数形式。本文中使用eig（）。e为计算得到的特征值，本文利用argsort函数将其从大到小排列，并选取前t个特征值。选取特征值的个数，即代表了降维后的维度。选择主成分个数，也就是k的个数可以通过计算方差保留比例。本实验中per1表示该百分比。代码如下。

[python] view plain copy

1. M=(dot(I1.T,I1))/(m-1)#计算协方差矩阵  
2. e,EV=eig(M)#计算特征值和特征向量，e为特征值，EV为特征向量  
3. e1=argsort(e)#计算e从小到大大索引值  
4. e1=e1[::-1]  
5. Q=EV[:,e1[0:t]]#将样本点投影到选取的特征向量上，选取几列特征向量意味着降为几维，这里降为2维  
6. percent=0  
7. for i in range(0,t):  
8.     percent=percent+e[i]  
9. per=percent/(sum(e))#计算差分百分比  
10. per1=real(per)*100  

（3） 原数据投影降维
将中心化的原数据与最大的k个特征值对应的特征向量点乘，得到降维后的新数据。I1（m*n）与Q（n*t）相乘，得到（m*t）数据，原来n维数据变为t维。再将原始数据投影到该子空间中。代码如下

[python] view plain copy

1. lowD = dot(I1, Q)#降维后的数据  
2. pj_train = dot(lowD.T, I1)#数据投影  

最后用matplotlib画图表现二维数据在PCA之后如何降为1维。实验全部代码如下。

[python] view plain copy

1. #-*- coding:utf-8 -*-  
2. from numpy import *  
3. from scipy.linalg import *  
4. from pylab import *  
5. import matplotlib.pyplot as plt  
6. def testeig(I,t):#提取前t个特征值  
7.     m,n = I.shape  
8.     mean_n = I.mean(axis=0)  
9.     I1 = I-mean_n  
10.     M = (dot(I1.T,I1))/(m-1)  
11.     e,EV = eig(M)  
12.     e1 = argsort(e)  
13.     e1 = e1[::-1]  
14.     Q = EV[:,e1[0:t]]  
15.     per=0  
16.     for i in range(t):  
17.         per+=e[e1[i]]  
18.     percent = per / sum(e)  
19.     percent1 = real(percent)*100  
20.     lowD = dot(I1, Q)#降维后的数据  
21.     pj_train = dot(lowD.T, I1)#数据投影  
22.     return pj_train,percent1  
23.   
24. I=rand(2,10)  
25. pj_train,percent1=testeig(I,2)  
26. plot(I[0,:],I[1,:],'rs')  
27. X=[I[0,:],pj_train[0,:]]  
28. Y=[I[1,:],pj_train[1,:]]  
29. plot(pj_train[0,:],pj_train[1,:],'bo')  
30. plot(X,Y)  
31. title('Plotting:PCA')  
32. show()  
33. print percent1,'\n',pj_train  

红色为原始数据坐标，蓝色为降维后坐标。相同点用线相连接。