22机器学习开放基础课程--主成分分析和聚类

主成分分析和聚类

无监督学习是机器学习中的一类重要算法。与监督学习算法相比，无监督学习算法不需要对输入数据进行标记，即不需要给出标签或类别。此外，无监督学习算法还可以在没有辅助的情况下能够学习数据的内在关系。
由于数据不需要手动标记标签，这可以使许多无监督算法能够处理大量的数据，从而节省了大量的人力成本。但是，这也给其带来一个问题，因为无监督学习算法在学习时没有使用标签数据，所以一般情况下很难直观地去评估无监督学习算法的质量。

主成分分析 PCA

在无监督学习中最常见的任务之一是降维，也就是减少输入数据的维数。 为什么要降维呢？主要有一下几个原因：首先，降维可能有助于数据可视化，因为人是无法理解高维数据的，通常只能看明白二维或三维数据的可视化图像。其次，降维可以有效的解决维度灾难的问题，改善模型的训练效果。此外，降维还可以进行数据压缩，提升模型训练效率。
目前，主成分分析是最简单，最直观，最常用的降维方法之一，其工作原理主要是将的数据投影到一个正交特征子空间中。

from mpl_toolkits.mplot3d import Axes3D
from sklearn import datasets
from sklearn import decomposition
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')
sns.set(style='white')
%matplotlib inline
%config InlineBackend.figure_format = 'retina'

通过 scikit-learn 提供的数据集接口导入莺尾花数据集。

iris = datasets.load_iris()
X = iris.data
y = iris.target

为了直观地查看数据的分布，使用三维图画出莺尾花的数据分布图。

fig = plt.figure(1, figsize=(6, 5))
plt.clf()
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=48, azim=134)

plt.cla()

for name, label in [('Setosa', 0), ('Versicolour', 1), ('Virginica', 2)]:
    ax.text3D(X[y == label, 0].mean(),
              X[y == label, 1].mean() + 1.5,
              X[y == label, 2].mean(), name,
              horizontalalignment='center',
              bbox=dict(alpha=.5, edgecolor='w', facecolor='w'))
# 改变标签的顺序，让其与数据匹配
y_clr = np.choose(y, [1, 2, 0]).astype(np.float)
ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=y_clr,
           cmap=plt.cm.nipy_spectral)

ax.w_xaxis.set_ticklabels([])
ax.w_yaxis.set_ticklabels([])
ax.w_zaxis.set_ticklabels([])

现在让我们看看使用 PCA 是怎么样提高模型的识别性能的。同样先导入实验所用到的模块。

from sklearn.tree import DecisionTreeClassifier  # 导入决策树模型、
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, roc_auc_score  # 识别准确率计算函数

莺尾花数据是一个相对容易区分的数据。因此，为了使实验结果对比明显。选用简单的决策树模型来对莺尾花数据进行分类。

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.3,
                                                    stratify=y,
                                                    random_state=42)

# 决策树的深度设置为 2
clf = DecisionTreeClassifier(max_depth=2, random_state=42)
clf.fit(X_train, y_train)
preds = clf.predict_proba(X_test)
print('Accuracy: {:.5f}'.format(accuracy_score(y_test,
                                               preds.argmax(axis=1))))

从上面的结果可知。在不对数据进行处理的情况下，使用决策树模型对莺尾花数据进行分类的准确率为 0.88889。
现在使用 PCA 将莺尾花数据的维度降低到 2 维，然后画出降维后的数据分布图。

# 使用从 sklearn 提供的 PCA 接口，降数据降到 2 维
pca = decomposition.PCA(n_components=2)
X_centered = X - X.mean(axis=0)
pca.fit(X_centered)
X_pca = pca.transform(X_centered)

# 可视化 PCA 降维后的结果
plt.plot(X_pca[y == 0, 0], X_pca[y == 0, 1], 'bo', label='Setosa')
plt.plot(X_pca[y == 1, 0], X_pca[y == 1, 1], 'go', label='Versicolour')
plt.plot(X_pca[y == 2, 0], X_pca[y == 2, 1], 'ro', label='Virginica')
plt.legend(loc=0)

同样的方法，将降维后的莺尾花数据输入到决策树模型中。

# 训练集合测试集同时使用 PCA 进行降维
X_train, X_test, y_train, y_test = train_test_split(X_pca, y, test_size=.3,
                                                    stratify=y,
                                                    random_state=42)

clf = DecisionTreeClassifier(max_depth=2, random_state=42)
clf.fit(X_train, y_train)
preds = clf.predict_proba(X_test)
print('Accuracy: {:.5f}'.format(accuracy_score(y_test,
                                               preds.argmax(axis=1))))

从上面的结果可知，对数据进行 PCA 降维之后，决策树模型的识别准确率提升到了 0.91111。这说明了 PCA 确实可以有效改善大部分机器学习算法的准确性和计算效率。
那么降维后的每个主成分都来自于原始数据的哪些维度呢？让我们通过查看每个维度的方差百分比来解释这个问题。

for i, component in enumerate(pca.components_):
    print("{} component: {}% of initial variance".format(i + 1,
                                                         round(100 * pca.explained_variance_ratio_[i], 2)))
    print(" + ".join("%.3f x %s" % (value, name)
                     for value, name in zip(component,
                                            iris.feature_names)))

从上面结果可知，我们将 4 维的数据降为了 2 维数据。在降维后的数据中，第一维（也就是第一个成分）主要由原始数据的0.361∗ sepallength −0.085∗ sepalwidth +0.857∗ petallength +0.358∗ petalwidth 组成。
手写数字数据集
在上面的例子中，莺尾花的原始数据只有 4 个维度。为了验证 PCA 在其他高维数据同样可行。接下来，使用之前实验所接触到的手写数字体数据集再完成一个示例练习。先导入手写数字数据集。

digits = datasets.load_digits()
X = digits.data
y = digits.target

数据集中的每个手写数字都是由 8×8 矩阵表示，每个像素值的表示颜色强度。获取数据集的前 10 个数字。并对这 10 个手写数字体数据进行可视化。

plt.figure(figsize=(14, 5))
for i in range(10):
    plt.subplot(2, 5, i + 1)
    plt.imshow(X[i, :].reshape([8, 8]), cmap='gray')

手写数字体数据集的维度是 8×8 维的，即 64 维。只有将其降维减少到 2 维，才能对其进行可视化，以便查看其分布状态。

pca = decomposition.PCA(n_components=2)
X_reduced = pca.fit_transform(X)

print('Projecting %d-dimensional data to 2D' % X.shape[1])

plt.figure(figsize=(6, 5))
plt.scatter(X_reduced[:, 0], X_reduced[:, 1], c=y,
            edgecolor='none', alpha=0.7, s=40,
            cmap=plt.cm.get_cmap('nipy_spectral', 10))
plt.colorbar()
plt.title('MNIST. PCA projection')

在上图中，总共包含十种颜色，每一种颜色对应一个类别标签。
除了 PCA 之外，t-SNE 也是一种常用的降维算法。相比于 PCA， t-SNE 不具有线性约束。下面使用 t-SNE 来对手写数字体数据进行降维，并对降维后的数据进行可视化。

from sklearn.manifold import TSNE  # 导入 t-SNE

tsne = TSNE(random_state=17)
X_tsne = tsne.fit_transform(X)

plt.figure(figsize=(6, 5))
plt.scatter(X_tsne[:, 0], X_tsne[:, 1], c=y,
            edgecolor='none', alpha=0.7, s=40,
            cmap=plt.cm.get_cmap('nipy_spectral', 10))
plt.colorbar()
plt.title('MNIST. t-SNE projection')

从上面的图中可以看出，t-SNE 相比于 PCA，降维效果要好很多。 但是 t-SNE 也有一个缺点，就是其运行时间要大大超过 PCA。在实际使用中，可以根据具体任务需求来选择。
在对手写数据集进行降维时，如果要保留原始数据的 90％ 散度，应该将数据降到多少维呢？先来画出主成分与其所保留的原始数据散度的关系。

pca = decomposition.PCA().fit(X)

plt.figure(figsize=(6, 4))
plt.plot(np.cumsum(pca.explained_variance_ratio_), color='k', lw=2)
plt.xlabel('Number of components')
plt.ylabel('Total explained variance')
plt.xlim(0, 63)
plt.yticks(np.arange(0, 1.1, 0.1))
plt.axvline(21, c='b')
plt.axhline(0.9, c='r')
plt.show()

从上图中可以看出，要想保留原始数据 90% 的信息。需要保留 21 个主要组成成分。因此，需要将维度从 64 个特征减少到 21 个。

聚类

聚类是无监督学习算法中的一种。聚类背后的主要思想相对简单，往往可以根据数据点之间的距离来将同类样本聚合在一起。

为了更好的理解 K-Means 算法的原理，这里通过一个例子来进行说明。先构建出一个数据集并画图它的分布图，该数据集含有三个簇。

# 构造可分为三个簇的数据
X = np.zeros((150, 2))

np.random.seed(seed=42)
X[:50, 0] = np.random.normal(loc=0.0, scale=.3, size=50)
X[:50, 1] = np.random.normal(loc=0.0, scale=.3, size=50)

X[50:100, 0] = np.random.normal(loc=2.0, scale=.5, size=50)
X[50:100, 1] = np.random.normal(loc=-1.0, scale=.2, size=50)

X[100:150, 0] = np.random.normal(loc=-1.0, scale=.2, size=50)
X[100:150, 1] = np.random.normal(loc=2.0, scale=.5, size=50)

plt.figure(figsize=(5, 5))
plt.plot(X[:, 0], X[:, 1], 'bo')

开始动手实现 K-Means 算法， K-Means 算法的实现非常简单。按照上述的算法步骤逐步实现即可。在这里，我们使用欧几里德距离来衡量两个数据点之间的距离。当然，你也可以使用其他距离度量方式。

# 调用 Scipy 库的距离计算函数，
# 用于计算数据点之间的距离
from scipy.spatial.distance import cdist

# 随机初始化三个中心点
np.random.seed(seed=42)
centroids = np.random.normal(loc=0.0, scale=1., size=6)
centroids = centroids.reshape((3, 2))

cent_history = []
cent_history.append(centroids)

for i in range(3):
    # 计算每个点到中心的距离
    distances = cdist(X, centroids)
    # 获取数据别分到哪个簇
    labels = distances.argmin(axis=1)
    # 根据数据到每个簇质心的距离，标记这些点的类别
    centroids = centroids.copy()
    centroids[0, :] = np.mean(X[labels == 0, :], axis=0)
    centroids[1, :] = np.mean(X[labels == 1, :], axis=0)
    centroids[2, :] = np.mean(X[labels == 2, :], axis=0)

    cent_history.append(centroids)

可视化出 K-Means 算法的运行过程，以便更好地理解。

# 可视化 K 均值聚类步骤
plt.figure(figsize=(8, 8))
for i in range(4):
    distances = cdist(X, cent_history[i])
    labels = distances.argmin(axis=1)

    plt.subplot(2, 2, i + 1)
    plt.plot(X[labels == 0, 0], X[labels == 0, 1], 'bo', label='cluster #1')
    plt.plot(X[labels == 1, 0], X[labels == 1, 1], 'co', label='cluster #2')
    plt.plot(X[labels == 2, 0], X[labels == 2, 1], 'mo', label='cluster #3')
    plt.plot(cent_history[i][:, 0], cent_history[i][:, 1], 'rX')
    plt.legend(loc=0)
    plt.title('Step {:}'.format(i + 1))

从上图中可以看出，仅用两步模型就收敛了，速度非常之快。但该算法有一个缺点，就是它对聚类质心的初始位置的很敏感，也即是随机选取的初始质心值。但是，可以多次运行算法，然后平均所有质心结果。
K-均值算法中 K 值的选择
相比于分类和回归等监督学习任务，聚类算法没有使用含标签的数据。因此，聚类需要更有效的模型评价方法。 通常，当使用 K-Means 算法时，需要最小化观测数据点与其所在的簇的质心之间的平方距离之和。该值越小说明该类聚得越好。
希望数据点尽可能接近它们所属簇的质心。但也存在一个问题，当质心数（也就是 K 的值）等于样本数时，公式得到最优解，即J(C) 达到最小。此时的每个样本都会单独把自己作为一类。这显然是没有意义的。因为这结果相当于没有进行聚类。为了避免这种情况，我们应该定义一个函数，使得J(C) 下降得不那么快。使用公式描述如下：

from sklearn.cluster import KMeans  # 导入 K-均值聚类模型

求出 K 值得选择与J(Ck) 的关系，并画出它们的关系图。

inertia = []
for k in range(1, 8):
    kmeans = KMeans(n_clusters=k, random_state=1).fit(X)
    inertia.append(np.sqrt(kmeans.inertia_))
plt.plot(range(1, 8), inertia, marker='s')
plt.xlabel('$k$')
plt.ylabel('$J(C_k)$')

从上图中，可以看到当k 小于 3 时，J(Ck)下降得非常快。之后相对平稳。这意味着选择k 等于 3 最为合适。
K-均值存在的问题

在每次迭代中，需要度量两个簇之间的距离。一般可以使用以下几种方法：

第三个平均连接是计算时间最有效的，因为它不需要在每次合并后，重新计算距离。
凝聚聚类的聚类结果可以用树状图进行表示，以帮助识别算法应该在哪一步停止以获得最佳的聚类结果。

现在用凝聚聚类实现一下前面使用 K-Means 实现的聚类例子

from scipy.cluster import hierarchy
from scipy.spatial.distance import pdist

X = np.zeros((150, 2))
# 构建数据集
np.random.seed(seed=42)
X[:50, 0] = np.random.normal(loc=0.0, scale=.3, size=50)
X[:50, 1] = np.random.normal(loc=0.0, scale=.3, size=50)

X[50:100, 0] = np.random.normal(loc=2.0, scale=.5, size=50)
X[50:100, 1] = np.random.normal(loc=-1.0, scale=.2, size=50)

X[100:150, 0] = np.random.normal(loc=-1.0, scale=.2, size=50)
X[100:150, 1] = np.random.normal(loc=2.0, scale=.5, size=50)

# pdist将计算成对距离矩阵的上三角形
distance_mat = pdist(X)
# 连接 - 是一种凝聚算法的实现
Z = hierarchy.linkage(distance_mat, 'single')
plt.figure(figsize=(10, 5))
dn = hierarchy.dendrogram(Z, color_threshold=0.5)

上图显示的是一个凝聚聚类的树状图。从上图可以清晰的看到整个聚类的过程：先把每个样本作为一个簇，经过多次迭代最终聚为一个簇的过程。

聚类模型评价

聚类算法与许多的监督学习算法不同，因为聚类在训练时不需要数据标签，所以不能简单地使用监督学习的评价方法来进行评价。
聚类算法通常有内部和外部的评价标准。 外部评价指标是使用有关已知真实划分的信息，也就是要知道数据的真实类别。而内部指标不需要知道数据的真实类别，仅根据聚类结果来进行评估。
调整兰德指数
兰德指数是一种常用的聚类评价方法。在该评价方法中，需要知道数据的真实标签，也需要数据集的聚类结果。 设n是样本中数据点对数的数量。 设a表示在真实标签与聚类结果中都是同类别的观测点对数。b表示在真实标签与聚类结果中都是不同类别的观测点对数。兰德指数可以表示为：

from sklearn import metrics
from sklearn import datasets
import pandas as pd
from sklearn.cluster import KMeans, AgglomerativeClustering, AffinityPropagation, SpectralClustering

加载数据，并构建不同的聚类算法。

data = datasets.load_digits()
X, y = data.data, data.target

algorithms = []
algorithms.append(KMeans(n_clusters=10, random_state=1))
algorithms.append(AffinityPropagation())
algorithms.append(SpectralClustering(n_clusters=10, random_state=1,
                                     affinity='nearest_neighbors'))
algorithms.append(AgglomerativeClustering(n_clusters=10))

使用不同评价指标对不同的聚类算法进行评价。

data = []
for algo in algorithms:
    algo.fit(X)
    data.append(({
        'ARI': metrics.adjusted_rand_score(y, algo.labels_),
        'AMI': metrics.adjusted_mutual_info_score(y, algo.labels_, average_method='arithmetic'),
        'Homogenity': metrics.homogeneity_score(y, algo.labels_),
        'Completeness': metrics.completeness_score(y, algo.labels_),
        'V-measure': metrics.v_measure_score(y, algo.labels_),
        'Silhouette': metrics.silhouette_score(X, algo.labels_)}))

results = pd.DataFrame(data=data, columns=['ARI', 'AMI', 'Homogenity', 'Completeness',
                                           'V-measure', 'Silhouette'],
                       index=['K-means', 'Affinity', 'Spectral', 'Agglomerative'])
results

非监督学习应用练习

我们将使用三星提供的 Human Activity Recognition 活动识别数据集 。这些数据来自三星 Galaxy S3 手机的加速度计和陀螺仪。这些活动类型包括：走路，站立，躺下，坐着或爬楼梯。
我们首先假装不了解活动的类型，并尝试纯粹基于特征对样本进行聚类。然后，我们将确定身体活动类型的问题解决为分类问题。
我们先导入本次挑战可能会用到的模块和函数。

from sklearn.svm import LinearSVC
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans, AgglomerativeClustering, SpectralClustering
from sklearn import metrics
from matplotlib import pyplot as plt
import os
import numpy as np
import pandas as pd
import seaborn as sns
from tqdm.notebook import tqdm
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
plt.style.use(['seaborn-darkgrid'])
plt.rcParams['figure.figsize'] = (12, 9)
plt.rcParams['font.family'] = 'DejaVu Sans'

RANDOM_STATE = 17

然后，加载并读取数据集。

!wget -nc "https://labfile.oss.aliyuncs.com/courses/1283/samsung_HAR.zip"
!unzip -o "samsung_HAR.zip"

PATH_TO_SAMSUNG_DATA = "./samsung_HAR"
X_train = np.loadtxt(os.path.join(PATH_TO_SAMSUNG_DATA, "samsung_train.txt"))
y_train = np.loadtxt(os.path.join(PATH_TO_SAMSUNG_DATA,
                                  "samsung_train_labels.txt")).astype(int)

X_test = np.loadtxt(os.path.join(PATH_TO_SAMSUNG_DATA, "samsung_test.txt"))
y_test = np.loadtxt(os.path.join(PATH_TO_SAMSUNG_DATA,
                                 "samsung_test_labels.txt")).astype(int)

X_train.shape, X_test.shape, y_train.shape, y_test.shape

聚类挑战中，我们不需要目标向量。所以，下面合并训练和测试样本：

X = np.vstack([X_train, X_test])
y = np.hstack([y_train, y_test])

X.shape, y.shape

我们尝试查看类别数据中唯一值有哪些，并得到数据集目标类别的数量 n_classes。

n_classes = np.unique(y).size
np.unique(y)

根据 数据集说明，我们可以得知类别数值表示的含义：
1 – walking
2 – walking upstairs
3 – walking downstairs
4 – sitting
5 – standing
6 – laying down
然后，尝试使用 StandardScaler 完成特征数据的规范化：

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
X_scaled.shape

接下来，我们将应用主成分分析 PCA 降维来缩减特征的数量。
如果使用 PCA 降维，并使得缩减后的数据保留原始数据 90% 的方差，那么得到的特征数量是多少个？

pca = PCA(n_components=0.9, random_state=RANDOM_STATE).fit(X_scaled)
X_pca = pca.transform(X_scaled)
X_pca.shape

降维后的第一主成分涵盖了多大比例的方差？舍入到最接近的百分比。

round(float(pca.explained_variance_ratio_[0] * 100))

请绘制前两个主成分特征的二维散点图，并使用数据已知类别进行着色。

plt.scatter(X_pca[:, 0], X_pca[:, 1], c=y, s=20, cmap='viridis')

接下来，请使用 KMeans 聚类方法对 PCA 降维后的数据进行聚类操作。于此同时，这里建议聚集为 6 类，实际上在正常情况下我们一般并不会知道聚为几类。
建立 KMeans 聚类，设置 n_clusters=n_classes, n_init=100 以及前面提供的 RANDOM_STATE。

kmeans = KMeans(n_clusters=n_classes, n_init=100,
                random_state=RANDOM_STATE, n_jobs=1)
kmeans.fit(X_pca)
cluster_labels = kmeans.labels_

同样，请使用 PCA 前面 2 个主成分特征绘制二维图像，但这一次使用 KMeans 聚类标签进行着色。

plt.scatter(X_pca[:, 0], X_pca[:, 1], c=cluster_labels, s=20, cmap='viridis')

对比查看原始标签分布和 KMeans 聚类标签的分布之间的不同之出。接下来，我们分别查看每个原始类别都被 KMeans 聚类划分成了那几个簇。

tab = pd.crosstab(y, cluster_labels, margins=True)
tab.index = ['walking', 'going up the stairs',
             'going down the stairs', 'sitting', 'standing', 'laying', 'all']
tab.columns = ['cluster' + str(i + 1) for i in range(6)] + ['all']
tab

表格行名称是原始标签，而列名称则对应了 KMeans 聚类后的簇序号。可以看到，几乎每一个原始类别都被重新聚为了几个分散簇。
这里，我们使用某一原始类别被 KMeans 聚类后的最大数量簇，除以原始类别总数来表征聚类的分散程度。得到的数值越小，即代表聚类后的簇越分散。
例如，所示的表格。walking 被重新划分为 cluster1, cluster4, cluster5。其中，数量最多的是 cluster1，则分散程度：

请计算各原始类别聚类后的分散程度，并按照分散程度由大到小排序。

pd.Series(tab.iloc[:-1, :-1].max(axis=1).values /
          tab.iloc[:-1, -1].values, index=tab.index[:-1]).sort_values()

接下来，通过求解观测数据点与其所在的簇的质心之间的平方距离之和来选择本次数据的最佳聚类 K 值。
选择 KMeans 聚类合适的聚类 K 值。

inertia = []
for k in tqdm_notebook(range(1, n_classes + 1)):
    kmeans = KMeans(n_clusters=k, n_init=100,
                    random_state=RANDOM_STATE, n_jobs=1).fit(X_pca)
    inertia.append(np.sqrt(kmeans.inertia_))

plt.plot(range(1, 7), inertia, marker='s')

d = {}
for k in range(2, 6):
    i = k - 1
    d[k] = (inertia[i] - inertia[i + 1]) / (inertia[i - 1] - inertia[i])
d

不出意外时，这里的最佳聚类 K 值应该为 2，也就是将数据聚集为 2 类。
接下来，尝试使用 Agglomerative 聚类，并计算二者（KMeans）的兰德指数 ARI 值。
这里使用 n_clusters=n_classes 时得到的 KMeans 聚类结果即可。

ag = AgglomerativeClustering(n_clusters=n_classes, linkage='ward').fit(X_pca)
print('KMeans: ARI =', metrics.adjusted_rand_score(y, cluster_labels))
print('Agglomerative CLustering: ARI =',
      metrics.adjusted_rand_score(y, ag.labels_))

一般情况下，ARI 值越高代表聚类效果越好。