python 实现PCA

Principal Component Analysis

PCA 一种非监督学习，常被用于数据降维、有损数据压缩、特征抽取、数据可视化。它也被称为Karhunen-Loève变换。

PCA的思想是将n维特征映射到k维空间上。k<n，这k维特征是全新的正交特征，是重新构造出来的k维特征，而不是简单地从n维特征中去除其余n−k维特征。

PCA算法流程

算法输入：数据集X-m*n

• 按列计算数据集X的均值X-mean，然后令Xnew=X − Xmean；
• 求解矩阵Xnew的协方差矩阵，并将其记为Cov；
• 计算协方差矩阵Cov的特征值和相应的特征向量；
• 将特征值按照从大到小的排序，选择其中最大的k个，然后将其对应的k个特征向量分别作为列向量组成特征向量矩阵W-n*k;
• 计算XnewW，即将数据集Xnew投影到选取的特征向量上，这样就得到了我们需要的已经降维的数据集XnewW。

注意，计算一个n*n矩阵的完整的特征向量分解的时间复杂度为 O(n3) 。如果我们将数据集投影到前 k 个主成分中，那么我们只需寻找前 kk个特征值和特征向量。这可以使用更高效的方法得到,例如幂方法(power method) (Golub and Van Loan, 1996)，它的时间复杂度为 O(kn2)，或者我们也可以使用 EM 算法。

PCA优缺点

优点：

• 它是无监督学习，完全无参数限制的。在PCA的计算过程中完全不需要人为的设定参数或是根据任何经验模型对计算进行干预，最后的结果只与数据相关，与用户是独立的。

• 用PCA技术可以对数据进行降维，同时对新求出的“主元”向量的重要性进行排序，根据需要取前面最重要的部分，将后面的维数省去，可以达到降维从而简化模型或是对数据进行压缩的效果。同时最大程度的保持了原有数据的信息。

• 各主成分之间正交，可消除原始数据成分间的相互影响。

• 计算方法简单，易于在计算机上实现。

缺点：

• 如果用户对观测对象有一定的先验知识，掌握了数据的一些特征，却无法通过参数化等方法对处理过程进行干预，可能会得不到预期的效果，效率也不高。
• 贡献率小的主成分往往可能含有对样本差异的重要信息。
• 特征值矩阵的正交向量空间是否唯一有待讨论。
• 在非高斯分布的情况下，PCA方法得出的主元可能并不是最优的，此时在寻找主元时不能将方差作为衡量重要性的标准。

Python 调用numpy 实现PCA

import numpy as np
import matplotlib.pyplot as plt

x=np.array([2.5,0.5,2.2,1.9,3.1,2.3,2,1,1.5,1.1])
y=np.array([2.4,0.7,2.9,2.2,3,2.7,1.6,1.1,1.6,0.9])
# s1 求平均值以及标准化 normalization
mean_x=np.mean(x)
mean_y=np.mean(y)
# 减去均值 为了数据特征标准化 
# 在图像处理中，常对图像 减去他的平均亮度，图像的整体明亮程度并不会影响图像中存在的是什么物体
scaled_x=x-mean_x
scaled_y=x-mean_y
#构造矩阵
data=np.matrix([[scaled_x[i],scaled_y[i]] for i in range(len(scaled_x))])

plt.plot(scaled_x,scaled_y,'o')    


# s2 求协方差矩阵(Covariance Matrix)
cov=np.cov(scaled_x,scaled_y)
print(cov,end="\n")       

# 除了求协方差矩阵，还能求散度矩阵Scatter Matrix
# 由于我们之前已经做过normalization，因此对于我们来说，
# 这个矩阵就是 data*data的转置矩阵。
scatter_matrix=np.dot(np.transpose(data),data)
print(scatter_matrix,end="\n")

# 其实协方差矩阵和散度矩阵关系密切，散度矩阵 就是协方差矩阵乘以（总数据量-1）
# 因此他们的特征根和特征向量是一样的。这里值得注意的一点就是，
# 散度矩阵是SVD奇异值分解的一步，因此PCA和SVD是有很大联系的

# s3 求协方差矩阵的特征根和特征向量
# 用numpy计算特征根eig_val和特征向量eig_vec
eig_val, eig_vec = np.linalg.eig(cov)
print("特征值")
print(eig_val,end="\n")
print("特征向量")
print(eig_vec,end="\n")

xmin ,xmax = scaled_x.min(), scaled_x.max()
ymin, ymax = scaled_y.min(), scaled_y.max()

dx = (xmax - xmin) * 0.2
dy = (ymax - ymin) * 0.2
plt.xlim(xmin - dx, xmax + dx)
plt.ylim(ymin - dy, ymax + dy)

# 其中红线、蓝线就是特征向量
plt.plot([eig_vec[:,0][0],0],[eig_vec[:,0][1],0],color='red')
plt.plot([eig_vec[:,1][0],0],[eig_vec[:,1][1],0],color='blue')

# 注意：
# 1.特征向量之间是正交的，PCA其实就是利用特征向量的这个特点，重新构建新的空间体系
# 2.特征向量代表着数据的pattern(模式),比如一条代表着y随着x的增大而增大的趋势，
# 而另外一条，则是代表数据也有该方面的变化。所以特征向量的命名是很科学的，他代表着矩阵的特征。

# 如果我们将数据直接乘以特征向量矩阵，那么其实我们只是以特征向量为基底，重新构建了空间
new_data=np.transpose(np.dot(eig_vec,np.transpose(data)))
print("将数据投影到特征向量空间之后的值")
print(new_data)
plt.plot(new_data[:,0],new_data[:,1],'^',color='blue')
# s4 得到特征值和特征向量之后，我们可以根据特征值的大小，从大到小的选择K个特征值对应的特征向量。
eig_pairs = [(np.abs(eig_val[i]), eig_vec[:,i]) for i in range(len(eig_val))]
print("排序前")
print(eig_pairs)
eig_pairs.sort(reverse=True)
print("排序后")
print(eig_pairs)

# 从eig_pairs选取前k个特征向量就行。这里，我们只有两个特征向量，选一个最大的。
feature=eig_pairs[0][1]

# s5 转化得到降维的数据
# 将原来的数据乘以经过筛选的特征向量组成的特征矩阵之后，就可以得到新的数据
print("new data:")
new_data_reduced=np.transpose(np.dot(feature,np.transpose(data)))
print(new_data_reduced)
plt.plot(new_data_reduced[:,0],[1.2]*10,'*',color='green')
plt.show()
# 绿色的五角星是PCA处理过后得到的一维数据，为了能跟以前的图对比，将他们的高度定位1.2，
# 其实就是红色圆点投影到蓝色线之后形成的点。

# 这就是PCA,通过选择特征根向量，形成新的坐标系，然后数据投影到这个新的坐标系，
# 在尽可能少的丢失信息的基础上实现降维