主成分分析是我十分喜欢的一个降维方法，在之前学习并总结过其matlab的实现方法，详情见傻瓜攻略（一）——MATLAB主成分分析代码及结果分析实例一文，最近学习自编码器的相关知识，经常看见的一句话是“自编码器可以用于特征降维，类似主成分分析PCA，但是其相比PCA其性能更强”，作为一个历史较长的方法，被新方法超越是十分常见的，但是不可否认PCA还是具有其独特的功能，在许多研究领域被广泛应用。
本文的基本内容来自吴茂贵《Python深度学习：基于TensorFlow》3.9小节，用Python实现鸢尾花数据的主成分分析，本文尝试理解代码并改进，可以对自己的数据进行分析。
主成分分析原理
主成分分析法(PCA)是一种高效处理多维数据的多元统计分析方法，将主成分分析用于多指标（变量）的综合评价较为普遍。笔者自从本科学习数学建模就开始接触该方法，但是一直没有系统地整理过，借这个机会总结一下，以备不时之需。
该方法的基本思想是运用较少的变量去解释原始数据中的大部分变异，通过对原始数据相关矩阵内部结构关系的分析和计算，产生一系列互不相关的新变量。根据需要从中选取比原始变量个数少的几个新变量，这些新的变量就是所谓的主成分，它们能够充分解释原始数据的变化。因此，主成分分析法本质上是一种降维方法，也多被用于高维数据的降维处理。
基本步骤
针对待分析数据X（n个特征），PCA的具体步骤如下：
（1）对向量X进行去中心化；
（2）计算X的协方差矩阵，自由度选择0或1；
（3）计算协方差矩阵的特征值和特征向量；
（4）选取最大的k个特征值及其特征向量（k<n）；
（5）用X与选取的特征向量相乘，得到k个主成分，即降维结果。
通过以上步骤，将原数据的n个特征减少到k个。
python代码
## pca特征降维
# 导入相关模块
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from numpy.linalg import eig
from sklearn.datasets import load_iris
# 导入数据
iris = load_iris() # 150*4的矩阵，4个特征，行是个数，列是特征
X = iris.data
k = 2 #选取贡献最大的前2个特征作为主成分，根据实际情况灵活选取
# Standardize by remove average通过去除平均值进行标准化
X = X - X.mean(axis=0)
# Calculate covariance matrix:计算协方差矩阵：
X_cov = np.cov(X.T, ddof=0)
# Calculate eigenvalues and eigenvectors of covariance matrix
# 计算协方差矩阵的特征值和特征向量
eigenvalues, eigenvectors = eig(X_cov)
# top k large eigenvectors选取前k个特征向量
klarge_index = eigenvalues.argsort()[-k:][::-1]
k_eigenvectors = eigenvectors[klarge_index]
# X和k个特征向量进行点乘
X_pca = np.dot(X, k_eigenvectors.T)
运行结果
相比于matlab代码输出的结果，这里的直观结果较少，可以通过print函数显示相关计算结果。下面按照主成分分析的步骤进行说明。
首先是数据标准化，主成分分析的结果直接受指标量纲的影响。由于各指标的单位可能不一样，因此进行量化评分得到的数据值大小也是不同的。如埋深的单位是米，相应指标在0.8到2.0之间，而人口密度指标的数据值在280左右，这样会导致分析结果的不准确。因此数据的标准化是主成分分析的前提条件，所以实际中可以先把各指标的数据标准化。这里通过去除平均值进行标准化。
完成数据的标准化后，对所得结果计算得到标准化数据的协方差矩阵（标准化数据的协方差矩阵即相关系数矩阵，保存在X_cov中）。相关系数代表了不同指标之间的相关程度，绝对值越大代表相关性越高。相关性较高的变量之间存在信息上的重叠，信息重叠在很大程度上会影响评价结果的客观性，因此相关性矩阵可以证明进行主成分分析的必要性。
由相关系数矩阵可以计算出特征值与特征向量，计算得到与指标数量n相等的n个待选主成分。n个特征值代表了n个待选主成分对最终评价结果的贡献程度，特征值保存在eigenvalues中，从大到小排列。待选主成分的特征向量为n*n的矩阵，表示待选主成分和相应的原始数据的相关关系，其绝对值越大，则待选主成分对该指标的代表性越大。特征向量保存在eigenvectors中，每一列代表一个特征向量，对应一个主成分。
在主成分的选取上，对应的特征值大小是一个重要衡量因素，普遍的做法是保存特征值要大于1的主成分，舍弃特征值小于1 的主成分，因此最终的主成分个数会小于指标个数n。也可以根据贡献度大小，累计贡献度达到某个程度，不同标准有70%以上，85%以上或其他。这里选取所有特征值大于1的主成分，选取的主成分个数为k，k个主成分的特征向量保存在k_eigenvectors中。
原始数据X与选取的特征向量相乘，得到k个主成分，即降维结果，即pca函数的调用结果，保存在X_pca中。
个人数据的降维
个人数据的形式一般为矩阵，如下是18个输油管段在10个指标上的表现，即一个18*10的矩阵。将其保存到空白txt文件中并保存，作为程序的原始数据。下面的数据来源为论文《基于主成分-聚类分析法的管道风险评价方法》。
1.5 7.1 280.0 424.0 4.0 430.0 90.0 112.0 459.0 453.0
1.5 7.1 280.0 424.0 4.0 431.0 92.0 112.0 459.0 453.0
1.5 8.7 280.0 424.0 4.0 431.0 92.0 112.0 459.0 453.0
1.5 8.7 280.0 424.0 4.0 431.0 92.0 112.0 459.0 448.0
1.5 8.7 280.0 424.0 4.0 431.0 90.0 112.0 455.0 453.0
2.0 10.3 279.0 425.0 2.0 431.0 90.0 110.0 458.0 449.0
2.0 10.3 277.0 425.0 2.0 431.0 92.0 112.0 458.0 449.0
0.8 10.3 277.0 425.0 2.0 431.0 92.0 112.0 454.0 448.0
2.0 10.3 277.0 425.0 2.0 431.0 92.0 112.0 454.0 448.0
2.0 10.3 277.0 425.0 2.0 431.0 92.0 112.0 458.0 448.0
2.0 10.3 277.0 425.0 2.0 431.0 92.0 377.0 456.0 448.0
2.0 10.3 279.0 425.0 2.0 431.0 92.0 377.0 454.0 449.0
2.0 10.3 279.0 425.0 2.0 431.0 90.0 377.0 456.0 446.0
2.0 10.3 279.0 425.0 2.0 431.0 90.0 112.0 454.0 449.0
2.0 10.3 279.0 425.0 4.0 431.0 92.0 112.0 458.0 449.0
2.0 10.3 279.0 425.0 4.0 430.0 92.0 112.0 458.0 449.0
2.0 10.3 274.0 425.0 4.0 430.0 92.0 112.0 458.0 449.0
1.6 10.3 279.0 424.0 4.0 430.0 92.0 112.0 458.0 449.0
将以上矩阵保存在pca.txt文档中，和.py文件放在同一文件夹中，将以上代码中的
iris = load_iris()
X = iris.data # 150*4的矩阵，4个特征，行是个数，列是特征
替换为
X =np.loadtxt('./pca.txt',dtype=np.float32)
就可以实现以上数据的降维了，还是得到2个新的特征，具体的个数需要根据特征值的贡献率来确定。
绘制贡献率图像
在原代码后面添加以下代码，就可以绘制各个特征值的贡献率，
## 绘制贡献率图像
tot=sum(eigenvalues) # 计算特征值的和
var_exp=[(i / tot) for i in sorted(eigenvalues,reverse=True)] # 按照降序排列特征值，并计算贡献率
cum_var_exp=np.cumsum(var_exp) #累计贡献度
plt.bar(range(1,len(var_exp)+1),var_exp,alpha=0.5,align='center',label='individual var') # 绘制柱状图，
plt.step(range(1,len(var_exp)+1),cum_var_exp,where='mid',label='cumulative var') # 绘制阶梯图
plt.ylabel('variance rtion') # 纵坐标
plt.xlabel('principal components') # 横坐标
plt.legend(loc='best') # 图例位置，右下角
plt.show()
鸢尾花数据的绘制结果如下
笔者的其他博客，欢迎大家阅读学习，共同进步
傻瓜攻略（一）——MATLAB主成分分析代码及结果分析实例
286
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
