举例学习视觉化的主成分分析

主成分分析（PCA）是一种无监督的机器学习技术。主成分分析最流行的用途也许是降维。除了将PCA作为一种数据准备技术，我们还可以用它来帮助数据的可视化。一张图片胜过千言万语。随着数据的可视化，我们更容易获得一些洞察力，并决定我们机器学习模型的下一步。

在本教程中，你将发现如何使用PCA将数据可视化，以及使用可视化来帮助确定降维的参数。

完成本教程后，你将知道。

• 如何使用可视化的高维数据
• 什么是PCA的解释方差
• 从高维数据的PCA结果中直观地观察解释方差

让我们开始吧。

教程概述

本教程分为两部分；它们是。

• 高维数据的散点图
• 解释方差的可视化

前提条件

对于本教程，我们假设你已经熟悉了。

• 如何在Python中从头开始计算主成分分析（PCA）
  
• Python中用于降维的主成分分析

高维数据的散点图

可视化是从数据中获得洞察力的一个关键步骤。我们可以从可视化中了解到是否可以观察到某种模式，从而估计出哪种机器学习模型是合适的。

描绘二维的事物是很容易的。通常情况下，带有X轴和Y轴的散点图是二维的。将事物描绘成三维是有点挑战，但并非不可能。例如，在matplotlib中，可以绘制三维图。唯一的问题是在纸上或屏幕上，我们需要一次只能在一个视口或投影上看一个三维图。在matplotlib中，这是由仰角和方位角的程度控制的。描绘四维或五维的事物是不可能的，因为我们生活在一个三维的世界里，不知道这么高维度的事物会是什么样子。

这就是像PCA这样的降维技术发挥作用的地方。我们可以将维度降低到2或3，这样我们就可以将其可视化。让我们从一个例子开始。

我们从葡萄酒数据集开始，它是一个有13个特征和3个类别的分类数据集。有178个样本。

from sklearn.datasets import load_wine
winedata = load_wine()
X, y = winedata['data'], winedata['target']
print(X.shape)
print(y.shape)

(178, 13)
(178,)

在这13个特征中，我们可以任意挑选两个，然后用matplotlib绘图（我们用c 参数对不同的类进行着色）。

...
import matplotlib.pyplot as plt
plt.scatter(X[:,1], X[:,2], c=y)
plt.show()

或者我们也可以选择任何三个，然后用3D显示。

...
ax = fig.add_subplot(projection='3d')
ax.scatter(X[:,1], X[:,2], X[:,3], c=y)
plt.show()

但这些并不能揭示数据的样子，因为大部分的特征都没有显示出来。我们现在求助于主成分分析。

...
from sklearn.decomposition import PCA
pca = PCA()
Xt = pca.fit_transform(X)
plot = plt.scatter(Xt[:,0], Xt[:,1], c=y)
plt.legend(handles=plot.legend_elements()[0], labels=list(winedata['target_names']))
plt.show()

在这里，我们通过PCA将输入数据X ，转化为Xt 。我们只考虑前两列，其中包含最多的信息，并将其绘制成二维图。我们可以看到，紫色类是相当有特色的，但仍有一些重叠的地方。但是如果我们在PCA之前对数据进行缩放，结果就会不同。

...
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
pca = PCA()
pipe = Pipeline([('scaler', StandardScaler()), ('pca', pca)])
Xt = pipe.fit_transform(X)
plot = plt.scatter(Xt[:,0], Xt[:,1], c=y)
plt.legend(handles=plot.legend_elements()[0], labels=list(winedata['target_names']))
plt.show()

因为PCA对尺度很敏感，如果我们把每个特征归一化，StandardScaler ，我们可以看到更好的结果。在这里，不同的类别更加鲜明。通过观察这个图，我们确信像SVM这样的简单模型可以对这个数据集进行高精确度的分类。

把这些放在一起，下面是生成可视化的完整代码。

from sklearn.datasets import load_wine
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import matplotlib.pyplot as plt

# Load dataset
winedata = load_wine()
X, y = winedata['data'], winedata['target']
print("X shape:", X.shape)
print("y shape:", y.shape)

# Show any two features
plt.figure(figsize=(8,6))
plt.scatter(X[:,1], X[:,2], c=y)
plt.xlabel(winedata["feature_names"][1])
plt.ylabel(winedata["feature_names"][2])
plt.title("Two particular features of the wine dataset")
plt.show()

# Show any three features
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(projection='3d')
ax.scatter(X[:,1], X[:,2], X[:,3], c=y)
ax.set_xlabel(winedata["feature_names"][1])
ax.set_ylabel(winedata["feature_names"][2])
ax.set_zlabel(winedata["feature_names"][3])
ax.set_title("Three particular features of the wine dataset")
plt.show()

# Show first two principal components without scaler
pca = PCA()
plt.figure(figsize=(8,6))
Xt = pca.fit_transform(X)
plot = plt.scatter(Xt[:,0], Xt[:,1], c=y)
plt.legend(handles=plot.legend_elements()[0], labels=list(winedata['target_names']))
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.title("First two principal components")
plt.show()

# Show first two principal components with scaler
pca = PCA()
pipe = Pipeline([('scaler', StandardScaler()), ('pca', pca)])
plt.figure(figsize=(8,6))
Xt = pipe.fit_transform(X)
plot = plt.scatter(Xt[:,0], Xt[:,1], c=y)
plt.legend(handles=plot.legend_elements()[0], labels=list(winedata['target_names']))
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.title("First two principal components after scaling")
plt.show()

如果我们在不同的数据集上应用同样的方法，如MINST手写数字，散点图没有显示出明显的边界，因此它需要一个更复杂的模型，如神经网络来进行分类。

from sklearn.datasets import load_digits
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
import matplotlib.pyplot as plt

digitsdata = load_digits()
X, y = digitsdata['data'], digitsdata['target']
pca = PCA()
pipe = Pipeline([('scaler', StandardScaler()), ('pca', pca)])
plt.figure(figsize=(8,6))
Xt = pipe.fit_transform(X)
plot = plt.scatter(Xt[:,0], Xt[:,1], c=y)
plt.legend(handles=plot.legend_elements()[0], labels=list(digitsdata['target_names']))
plt.show()

可视化解释方差

PCA在本质上是通过其线性组合重新排列特征。因此，它被称为一种特征提取技术。PCA的一个特点是，第一个主成分拥有关于数据集的最多信息。第二个主成分比第三个主成分的信息量更大，以此类推。

为了说明这个想法，我们可以分步从原始数据集中去除主成分，看看数据集是什么样子的。让我们考虑一个特征较少的数据集，并在一个图中显示两个特征。

from sklearn.datasets import load_iris
irisdata = load_iris()
X, y = irisdata['data'], irisdata['target']
plt.figure(figsize=(8,6))
plt.scatter(X[:,0], X[:,1], c=y)
plt.show()

这是虹膜数据集，它只有四个特征。这些特征的尺度是相当的，因此我们可以跳过标度器。对于4个特征的数据，PCA最多可以产生4个主成分。

...
pca = PCA().fit(X)
print(pca.components_)

[[ 0.36138659 -0.08452251  0.85667061  0.3582892 ]
 [ 0.65658877  0.73016143 -0.17337266 -0.07548102]
 [-0.58202985  0.59791083  0.07623608  0.54583143]
 [-0.31548719  0.3197231   0.47983899 -0.75365743]]

例如，第一行是第一个主轴，在此基础上产生了第一个主成分。对于具有$p=(a,b,c,d)$特征的任何数据点，由于主轴由向量$v=(0.36,-0.08,0.86,0.36)$表示，该数据点的第一主成分在主轴上的值为$0.36\\times a - 0.08\\times b + 0.86\\times c + 0.36\\times d$。使用向量点积，这个值可以用

因此，将数据集$X$作为150 $/times$ 4矩阵（150个数据点，每个有4个特征），我们可以通过矩阵-向量乘法将每个数据点映射到这个主轴上的值：

，结果就是一个长度为150的向量。现在，如果我们从每个数据点中去除沿主轴向量的相应值，那将是
$X - (X \\times v) \\times v^T$
其中转置的向量 $v^T$ 是一行，$X\\times v$ 是一列。乘积$(X\\times v)\\times v^T$遵循矩阵-矩阵乘法，结果是一个$150times 4$矩阵，与$X$的维度相同。

如果我们绘制$(X\\times v) \\times v^T$的前两个特征，看起来是这样的。

...
# Remove PC1
Xmean = X - X.mean(axis=0)
value = Xmean @ pca.components_[0]
pc1 = value.reshape(-1,1) @ pca.components_[0].reshape(1,-1)
Xremove = X - pc1
plt.scatter(Xremove[:,0], Xremove[:,1], c=y)
plt.show()

numpy数组Xmean ，就是把X 的特征移到以零为中心。这是PCA的要求。然后，数组value 是通过矩阵-向量乘法计算出来的。
数组value 是在主轴上映射的每个数据点的大小。因此，如果我们将这个值与主轴向量相乘，就会得到一个数组pc1 。从原始数据集X ，我们得到一个新的数组Xremove 。在图中我们观察到，散点图上的点溃散在一起，每个类别的聚类也没有以前那么明显。这意味着我们通过去除第一个主成分，删除了很多信息。如果我们再次重复同样的过程，这些点就会进一步崩溃。

...
# Remove PC2
value = Xmean @ pca.components_[1]
pc2 = value.reshape(-1,1) @ pca.components_[1].reshape(1,-1)
Xremove = Xremove - pc2
plt.scatter(Xremove[:,0], Xremove[:,1], c=y)
plt.show()

这看起来像一条直线，但实际上不是。如果我们再重复一次，所有的点都塌陷成一条直线。

...
# Remove PC3
value = Xmean @ pca.components_[2]
pc3 = value.reshape(-1,1) @ pca.components_[2].reshape(1,-1)
Xremove = Xremove - pc3
plt.scatter(Xremove[:,0], Xremove[:,1], c=y)
plt.show()

这些点都落在一条直线上，因为我们从只有四个特征的数据中删除了三个主成分。因此，我们的数据矩阵变成了等级1。你可以试着再重复一次这个过程，结果会是所有的点都塌陷成一个点。当我们删除主成分时，每一步所删除的信息量可以通过PCA的相应解释方差率来发现。

...
print(pca.explained_variance_ratio_)

[0.92461872 0.05306648 0.01710261 0.00521218]

这里我们可以看到，第一个成分解释了92.5%的方差，第二个成分解释了5.3%的方差。如果我们去掉前两个主成分，剩下的方差只有2.2%，因此从视觉上看，去掉两个成分后的图看起来像一条直线。事实上，当我们用上面的图检查时，不仅我们看到点是崩溃的，而且随着我们去除成分，X轴和Y轴上的范围也变小了。

在机器学习方面，我们可以考虑在这个数据集中只使用一个单一的特征进行分类，即第一主成分。我们应该期望达到与使用全套特征一样的不低于90%的原始准确性。

...
from sklearn.model_selection import train_test_split
from sklearn.metrics import f1_score
from collections import Counter

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33)
from sklearn.svm import SVC
clf = SVC(kernel="linear", gamma='auto').fit(X_train, y_train)
print("Using all features, accuracy: ", clf.score(X_test, y_test))
print("Using all features, F1: ", f1_score(y_test, clf.predict(X_test), average="macro"))

mean = X_train.mean(axis=0)
X_train2 = X_train - mean
X_train2 = (X_train2 @ pca.components_[0]).reshape(-1,1)
clf = SVC(kernel="linear", gamma='auto').fit(X_train2, y_train)
X_test2 = X_test - mean
X_test2 = (X_test2 @ pca.components_[0]).reshape(-1,1)
print("Using PC1, accuracy: ", clf.score(X_test2, y_test))
print("Using PC1, F1: ", f1_score(y_test, clf.predict(X_test2), average="macro"))

Using all features, accuracy:  1.0
Using all features, F1:  1.0
Using PC1, accuracy:  0.96
Using PC1, F1:  0.9645191409897292

了解解释方差的另一个用途是在压缩上。鉴于第一主成分的解释方差很大，如果我们需要存储数据集，我们可以只存储第一主轴上的投影值（$X/times v$），以及主轴的向量$v$。然后，我们可以通过乘以它们来近似再现原始数据集：

这样，我们只需要为每个数据点存储一个值，而不是为四个特征存储四个值。如果我们将投影值存储在多个主轴上，并将多个主成分相加，那么近似值会更准确。

把这些放在一起，下面是生成可视化的完整代码。

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.decomposition import PCA
from sklearn.metrics import f1_score
from sklearn.svm import SVC
import matplotlib.pyplot as plt

# Load iris dataset
irisdata = load_iris()
X, y = irisdata['data'], irisdata['target']
plt.figure(figsize=(8,6))
plt.scatter(X[:,0], X[:,1], c=y)
plt.xlabel(irisdata["feature_names"][0])
plt.ylabel(irisdata["feature_names"][1])
plt.title("Two features from the iris dataset")
plt.show()

# Show the principal components
pca = PCA().fit(X)
print("Principal components:")
print(pca.components_)

# Remove PC1
Xmean = X - X.mean(axis=0)
value = Xmean @ pca.components_[0]
pc1 = value.reshape(-1,1) @ pca.components_[0].reshape(1,-1)
Xremove = X - pc1
plt.figure(figsize=(8,6))
plt.scatter(Xremove[:,0], Xremove[:,1], c=y)
plt.xlabel(irisdata["feature_names"][0])
plt.ylabel(irisdata["feature_names"][1])
plt.title("Two features from the iris dataset after removing PC1")
plt.show()

# Remove PC2
Xmean = X - X.mean(axis=0)
value = Xmean @ pca.components_[1]
pc2 = value.reshape(-1,1) @ pca.components_[1].reshape(1,-1)
Xremove = Xremove - pc2
plt.figure(figsize=(8,6))
plt.scatter(Xremove[:,0], Xremove[:,1], c=y)
plt.xlabel(irisdata["feature_names"][0])
plt.ylabel(irisdata["feature_names"][1])
plt.title("Two features from the iris dataset after removing PC1 and PC2")
plt.show()

# Remove PC3
Xmean = X - X.mean(axis=0)
value = Xmean @ pca.components_[2]
pc3 = value.reshape(-1,1) @ pca.components_[2].reshape(1,-1)
Xremove = Xremove - pc3
plt.figure(figsize=(8,6))
plt.scatter(Xremove[:,0], Xremove[:,1], c=y)
plt.xlabel(irisdata["feature_names"][0])
plt.ylabel(irisdata["feature_names"][1])
plt.title("Two features from the iris dataset after removing PC1 to PC3")
plt.show()

# Print the explained variance ratio
print("Explainedd variance ratios:")
print(pca.explained_variance_ratio_)

# Split data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33)

# Run classifer on all features
clf = SVC(kernel="linear", gamma='auto').fit(X_train, y_train)
print("Using all features, accuracy: ", clf.score(X_test, y_test))
print("Using all features, F1: ", f1_score(y_test, clf.predict(X_test), average="macro"))

# Run classifier on PC1
mean = X_train.mean(axis=0)
X_train2 = X_train - mean
X_train2 = (X_train2 @ pca.components_[0]).reshape(-1,1)
clf = SVC(kernel="linear", gamma='auto').fit(X_train2, y_train)
X_test2 = X_test - mean
X_test2 = (X_test2 @ pca.components_[0]).reshape(-1,1)
print("Using PC1, accuracy: ", clf.score(X_test2, y_test))
print("Using PC1, F1: ", f1_score(y_test, clf.predict(X_test2), average="macro"))

摘要

在本教程中，您发现了如何使用主成分分析实现数据的可视化。

具体而言，您学会了

• 使用PCA将高维数据集在二维上可视化
• 如何使用PCA维度中的绘图来帮助选择一个合适的机器学习模型
• 如何观察PCA的解释方差率
• 解释方差率对机器学习意味着什么