阅读 29 SEO

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

主成分分析和聚类

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

主成分分析 PCA

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

image.png

为了更好的理解 PCA 的工作原理,一起来看上面这幅图。试想一下,如果将坐标系中的数据点投影到 xx 轴(投影过后的数据即为 xx 轴上的一堆数据点),或者以同样的方法投影到 yy 轴,哪一个更容易区分呢?
很显然是 xx 轴。因为数据投影到 xx 轴之后相比于 yy 轴更分散,也就是数据样本的方差更大,从而也更容易区分。回到上图,如果要在坐标系中找到一根轴,使得投影到这根轴后的数据的方差更大,也就是说更分散,哪一根轴最好呢?显然,就是图中所标的椭圆轴方向。
PCA 的思想就是这样。如果将图中的数据都投影到椭圆轴上,则数据也就从二维降为了一维。
为了将数据的维数从n维降到k维,我们按照散度(也就是数据的分散程度)降低的顺序对轴列表进行排序,然后取出前k项。
现在开始计算原始数据 nn 维的散度值和协方差。根据协方差矩阵的定义,两个特征列的协方差矩阵计算公式如下:
image.png

上式中的μi表示第 i个特征的期望值。这里需要注意的是,协方差都是对称的,而且向量与其自身的协方差就等于其散度。当协方差为正时,说明X和Y是正相关关系;协方差为负时,说明X和Y是负相关关系;协方差为 0 时,说明X和Y是相互独立。
假定X 是观测矩阵,则协方差矩阵如下:
image.png

假定一个矩阵M的特征向量wi以及对应的特征向量为λi,则它们会满足下面公式:
image.png

一个样本X 的协方差矩阵可以表示为其转置矩阵X^T 与X 的乘积。根据 Rayleigh_quotient ,样本X 的最大方差位于协方差矩阵的最大特征值对应的特征向量上。也就是说想要保留一个矩阵的最大信息,我们只需要保留该矩阵的最大特征值所对应的特征向量所组成的矩阵即可,这个过程就是降维了。
因此,从数据中保留的主要成分就是与矩阵的顶部 kk 最大特征值对应的特征向量。
莺尾花数据集
上面主要讲述了主成分分析方法的原理,现在通过实例来加深理解。首先导入所有实验所用到的基本模块。

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 个。

聚类

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

image.png

K-Means 聚类
K-Means 算法是所有聚类算法中最流行且最简单的算法。 下面是它的工作原理:
选择认为最佳的类别数量k,即样本大概可以分为多少个簇。
在数据空间内随机初始化k 点为“质心”。
将每个观察数据点划分到于其最近的簇的质心的簇。
将质心更新为一个簇中所有数据点的中心。
重复步骤 3 和 4 步骤直到所有质心都相对稳定。

为了更好的理解 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) 下降得不那么快。使用公式描述如下:

image.png

为了更好的理解,先看看一个例子。因为 scikit-learn 也提供了各种聚类算法接口,而且使用这些接口有许多优点。例如:这些算法可以并行完成,有效减少了计算时间。所以在这里为了方便,直接使用 scikit-learn 提供的 K-Means 接口进行实验。

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-均值存在的问题

image.png

近邻传播算法
近邻传播也属于聚类算法中的一种。 与 K-Means 不同的是,这种方法不需要事先设置簇的数量,也就是 K 的值。 近邻传播算法的主要思想是根据观测点之间的相似性来对数据进行聚类。
image.png

image.png

谱聚类
谱聚类也是一种常用的聚类算法,其结合了上述的一些方法来创建一个性能更强的聚类方法。与其他聚类方法不同,谱聚类是一种图模型。谱聚类将观测点作为图模型的节点,将观测点之间的相似度看作为图模型中节点之间的边。在图模型中,如果两个观测点的边越短,则两个观测点就越相似。
在谱聚类中,可以通过定义观测点的相似性矩阵来存放数据样本中每两个观测点之间的相似性,通常也称为邻接矩阵。 这可以采用与近邻传播算法类似的方式来构建完成:
image.png

相似矩阵相当于描绘了一张图。当图模型构建完成时,就可以对其进行切割。切割原则就是较长的边给切开。换句话说就是把相似的观测点分为一类。 在形式上,这是一个 Normalized cuts 问题。
凝聚聚类
如果你不知道要将数据样本聚为多少个簇,可以考虑使用凝聚聚类算法。
凝聚聚类的算法流程很简单,如下:
首先将每个观测点都作为一个簇
然后按降序对每两个簇中心之间距离进行排序
取最近的两个相邻的簇并将它们合并为一个簇,然后重新计算簇中心
重复步骤 2 和 3 ,直到所有观测点都合并到一个簇中

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


image.png

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

现在用凝聚聚类实现一下前面使用 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表示在真实标签与聚类结果中都是不同类别的观测点对数。兰德指数可以表示为:

image.png

换句话说,兰德指数评估分割后的聚类结果和初始标签一致的比例。对于随机结果,RI 并不能保证分数接近 0 。为了实现在聚类结果是随机产生的情况下,兰德指数应该接近零,所以有必要缩放其大小,由此得到了调整兰德指数:
image.png

该评价方法是对称的,不依赖于标签排列顺序。 因此,该指数是样本不同划分之间距离的度量。 ARI 可取 [−1,1] 范围内的值。 其值越大表示聚类结果越好。
调整互信息
该评价指标与ARI 类似。 它也是对称的,不依赖于标签的具体值和排列顺序。 互信息MI 由 信息熵 函数定义,用来衡量真实数据分布与聚类结果分布的吻合程度。
同ARI 一样,这里也定义了调整互信息AMI 。 这可以避免随着聚类数量的增加而导致MI 指数增长。AMI 位于[0,1] 范围内。 其值越大意味着聚类结果越好。
同质性、完整性、V-mearure
从形式上来看,这些指标也是基于熵函数和条件熵函数来定义的,将聚类结果看作是离散分布:
image.png

其中K是聚类结果,C是原始数据。 因此,h评估的是每个簇是否由相同类别的数据点组成,c评估的是相同类别的数据点与所属簇的匹配程度。 这些指标不是对称的。两者都位于[0,1]范围内,而接近1的值表示更准确的聚类结果。 这些指标的值不像AMI 指标进行缩放过,因此取决于簇的数量。
当簇的数量很大且观测点的数量很小时,该评价方法也不会使评估值接近零。 在这种情况下,使用 ARI 可能会更合理。 但是,如果观测点超过 100 个,且簇的数量少于 10 个,此问题则可以忽略。
V−mearure 结合了h和c,为h和c的调和平均值。用来衡量两个划分结果的一致性。其公式如下:
image.png

轮廓系数
与上述评价指标相反,轮廓系数并不需要关于原始数据的真实标签信息。仅使用聚类结果来估计聚类的质量。设a是数据点与一个簇内其他观测点之间距离的平均值,b是观测点到最近簇的观测点的平均距离。则一个观测点的轮廓值为:
image.png

样本的轮廓值是所有样本轮廓值的平均值。该系数的取值范围为[−1,1]。轮廓系数值越高,表示聚类的结果越好。
因此,通常如果不知道数据可以分为几类,则可以通过获取最大化轮廓系数来确定最佳聚类数k。
一般情况下,都不需要我们去实现上述的评价方法,因为在 sklearn.metrics 接口中已经提供,只需导入即可。
最后,通过实验来看看这些指标是如何评价 MNIST 手写数字数据集的,先导入必要的模块:

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,则分散程度:


image.png

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

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 值越高代表聚类效果越好。

作者:Jachin111

原文链接:https://www.jianshu.com/p/5194fba43345

文章分类
后端
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容, 请发送邮件至 gxwowoo@163.com 举报,一经查实,本站将立刻删除。
相关推荐