基于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.stanford.edu/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进行计算，需先对其转置。其代码如下。

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

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

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

（3） 原数据投影降维

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

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

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

#-*- coding:utf-8 -*-
from numpy import *
from scipy.linalg import *
from pylab import *
import matplotlib.pyplot as plt
def testeig(I,t):#提取前t个特征值
m,n = I.shape
mean_n = I.mean(axis=0)
I1 = I-mean_n
M = (dot(I1.T,I1))/(m-1)
e,EV = eig(M)
e1 = argsort(e)
e1 = e1[::-1]
Q = EV[:,e1[0:t]]
per=0
for i in range(t):
per+=e[e1[i]]
percent = per / sum(e)
percent1 = real(percent)*100
lowD = dot(I1, Q)#降维后的数据
pj_train = dot(lowD.T, I1)#数据投影
return pj_train,percent1

I=rand(2,10)
pj_train,percent1=testeig(I,2)
plot(I[0,:],I[1,:],'rs')
X=[I[0,:],pj_train[0,:]]
Y=[I[1,:],pj_train[1,:]]
plot(pj_train[0,:],pj_train[1,:],'bo')
plot(X,Y)
title('Plotting:PCA')
show()
print percent1,'\n',pj_train

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