手把手教你使用PCA进行数据降维
对数据降维可以帮助我们提取数据集的主要信息，即将原始的高维特征空间压缩到低纬度的特征子空间。数据降维是用于提高计算效率的典型手段，另一个好处是也能够减小维度诅咒。
PCA(principal component analysis, 主成分分析)是一种被广泛使用的无监督的线性转换技术，主要用于降维。其他领域的应用还包括探索数据分析和股票交易的信号去噪，基因数据分析和基因表达。今天我们来结合代码学习一下PCA对数据降维的一个流程。
什么是PCA
PCA根据特征之间的相关性帮助我们确定数据中存在的模式。简而言之，PCA的目标是找到高维数据中最大方差的方向，并且将高维数据映射到一个新的子空间，这个子空间的方向不大于原始特征空间。新子空间的正交轴(主成分)可以被解释为原始空间的最大方差方向。下图
中是原始特征轴， 是主成分：
由于PCA方向极其容易受到数据中特征范围影响，所以在运用PCA前一定要做特征标准化，这样才能保证每维度特征的重要性等同。
PCA流程
1
数据标准化
我们使用一个叫做Wine数据集，先将原始Wine分割为训练集和测试集，然后标准化：
import pandas as pd
df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/'
'machine-learning-databases/wine/wine.data',
header=None)
df_wine.columns = ['Class label', 'Alcohol', 'Malic acid', 'Ash',
'Alcalinity of ash', 'Magnesium', 'Total phenols',
'Flavanoids', 'Nonflavanoid phenols', 'Proanthocyanins',
'Color intensity', 'Hue',
'OD280/OD315 of diluted wines', 'Proline']
df_wine.head()运行结果：
接下来是分割数据集
from distutils.version import LooseVersion as Version
from sklearn import __version__ as sklearn_version
if Version(sklearn_version) < '0.18':
from sklearn.cross_validation import train_test_split
else:
from sklearn.model_selection import train_test_split
X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
X_train, X_test, y_train, y_test = \
train_test_split(X, y, test_size=0.3, random_state=0)最后是标准化数据
from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.transform(X_test)上面的代码完成了PCA的第一步，对数据进行标准化。
2
构建协方差矩阵
我们再来看第二步：构建协方差矩阵。协方差矩阵是对称矩阵，d*d维度，其中d是原始数据的特征维度，协方差矩阵的每个元素是两两特征之间的协方差。
假设数据有三维度特征，协方差矩阵如下：
协方差矩阵的特征向量代表了主成分(最大方差的方向)，对应的特征值决定了特征向量绝对值的大小。在Wine数据集对应的13*13的协方差矩阵，我们会得到13个特征值。
NumPy提供了linalg.eig函数用于得到特征值和特征向量：
import numpy as np
cov_mat = np.cov(X_train_std.T)
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
print('\nEigenvalues \n%s' % eigen_vals)运行结果：
我们先使用np.cov方法得到数据的协方差矩阵，然后利用linalg.eig方法计算出特征向量(eigen_vecs)和特征值(eigen_vals)。
由于我们的目的是数据压缩，即降维，所以我们只将那些包含最多信息(方差)的特征向量(主成分)拿出来。什么样的特征向量包含的信息最多呢？这就要看特征值了，因为特征值定义了特征向量的大小，我们先对特征值进行降序排序，前k个特征值对应的特征向量就是我们要找
的主成分。
我们再学一个概念：方差解释率(variance explained ration)。一个特征值的方差解释率就是次特征值在特征值总和的占比。
利用NumPy提供的cumsum函数，我们可以计算累计解释方差和：
tot = sum(eigen_vals)
var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)import matplotlib.pyplot as plt
plt.bar(range(1, 14), var_exp, alpha=0.5, align='center',
label='individual explained variance')
plt.step(range(1, 14), cum_var_exp, where='mid',
label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.legend(loc='best')
plt.tight_layout()
# plt.savefig('./figures/pca1.png', dpi=300)
plt.show()
从上面的结果图我们可以看到第一个主成分占了近40%的方差(信息)，前两个主成分占了60%的方差。方差的物理含义是对值沿着特征轴的传播进行度量。
3
特征转换
在得到特征向量后，接下来我们就可以对原始特征进行转换了。本节我们先对特征值进行降序排序，然后用特征向量构建映射矩阵，最后用映射矩阵将原始数据映射到低维度特征子空间。
先对特征值排序：
# Make a list of (eigenvalue, eigenvector) tuples
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i])
for i in range(len(eigen_vals))]
# Sort the (eigenvalue, eigenvector) tuples from high to low
eigen_pairs.sort(key=lambda k: k[0], reverse=True)接下来，我们选择最大的两个特征值对应的特征向量，这里只用两个特征向量是为了下面画图方便，在实际运用PCA时，到底选择几个特征向量，要考虑到计算效率和分类器的表现两个方面(常用的选择方式是特征值子集要包含90%方差)：
w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
eigen_pairs[1][1][:, np.newaxis]))
print('Matrix W:\n', w)运行结果：
特征维度降到2维度后，我们就可以用散点图将数据可视化出来了：
X_train_pca = X_train_std.dot(w)
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']
for l, c, m in zip(np.unique(y_train), colors, markers):
plt.scatter(X_train_pca[y_train == l, 0],
X_train_pca[y_train == l, 1],
c=c, label=l, marker=m)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
# plt.savefig('./figures/pca2.png', dpi=300)
plt.show()
从上图可以看到，数据在x轴(第一主成分)上要比y轴(第二主成分)分布更广，这也符合方差解释率的结果。数据降维后，直觉上使用线性分类器就能够将数据分类。
scikit-learn中的PCA
上一小节我们详细讨论了PCA的步骤，在实际应用时，一般不会使用自己实现，而是直接调用sklearn中的PCA类，PCA类是另一个transformer类：我们先用训练集训练模型参数，然后统一应用于训练集和测试集。
下面我们就是用sklearn中的PCA类对Wine数据降维，然后调用逻辑回归模型分类，最后将决策界可视化出来：
from sklearn.decomposition import PCA
pca = PCA()
X_train_pca = pca.fit_transform(X_train_std)
pca.explained_variance_ratio_plt.bar(range(1, 14), pca.explained_variance_ratio_, alpha=0.5, align='center')
plt.step(range(1, 14), np.cumsum(pca.explained_variance_ratio_), where='mid')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.show()
pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.transform(X_test_std)plt.scatter(X_train_pca[:, 0], X_train_pca[:, 1])
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.show()
from matplotlib.colors import ListedColormap
def plot_decision_regions(X, y, classifier, resolution=0.02):
# setup marker generator and color map
markers = ('s', 'x', 'o', '^', 'v')
colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
cmap = ListedColormap(colors[:len(np.unique(y))])
# plot the decision surface
x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
np.arange(x2_min, x2_max, resolution))
Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
Z = Z.reshape(xx1.shape)
plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
plt.xlim(xx1.min(), xx1.max())
plt.ylim(xx2.min(), xx2.max())
# plot class samples
for idx, cl in enumerate(np.unique(y)):
plt.scatter(x=X[y == cl, 0],
y=X[y == cl, 1],
alpha=0.6,
c=cmap(idx),
edgecolor='black',
marker=markers[idx],
label=cl)from sklearn.linear_model import LogisticRegression
lr = LogisticRegression()
lr = lr.fit(X_train_pca, y_train)plot_decision_regions(X_train_pca, y_train, classifier=lr)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
# plt.savefig('./figures/pca3.png', dpi=300)
plt.show()执行上面的代码，我们应该得到如下的决策界：
如果你仔细观察我们自己实现的PCA得到的散点图和调用sklearn中PCA得到的散点图，两幅图看起来是镜像关系！这是由于NumPy和sklearn求解特征向量时计算的差异，如果你实在看不惯，只需要将其中一个得到的特征向量*(-1)即可。还要注意特征向量一般都要归一化。
参考：
Python Machine Learning Sebastian Raschka