AsanoSaki


夏天的海边有冰淇淋口味的海风~

KMeans聚类与PCA主成分分析

  • ~6.89K words
  1. 1. 二维数据K-Means聚类
  2. 2. 主成分分析PCA
  3. 3. 高维数据PCA后聚类
  4. 4. 高维数据聚类并计算与中心的相似度

介绍二维数据与高维数据的 K-Means 聚类算法以及高维数据的 PCA 主成分分析方法。

1. 二维数据K-Means聚类

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans, MiniBatchKMeans
from sklearn.decomposition import PCA


# 以center为中心产生随机分布的数据
def get_random_data(center, data_num, data_dim=2, fix_seed=True):
    if fix_seed:
        np.random.seed(0)
    # 产生-2~2的随机数
    offset = np.random.rand(data_num, data_dim) * 4 - 2
    return center + offset


centers = [[2, 8], [5, 2], [9, 6]]
center_num = len(centers)
data_num = 500  # 每个center产生的样本数量

data = np.zeros([data_num * center_num, 2])
for i, center in enumerate(centers):
    data[i * data_num:(i + 1) * data_num] = get_random_data(center, data_num)

# 显示数据
print(data.shape)  # (1500, 2)
plt.scatter(data[:, 0], data[:, 1], s=5, color='c')
plt.show()


# ---------- KMeans聚类 ----------
kms = KMeans(n_clusters=3, init='k-means++')  # 聚出3类
kms.fit(data)  # 模型拟合
centers = kms.cluster_centers_  # 计算聚类中心
labels = kms.labels_  # 聚类后每个样本的类别标签


# ---------- MiniBatchKMeans聚类 ----------
# mbk = MiniBatchKMeans(init='k-means++', n_clusters=3, batch_size=32, random_state=0)
# mbk.fit(data)
# centers = mbk.cluster_centers_
# labels = mbk.labels_

plt.figure(figsize=(8, 6))
plt.scatter(data[:, 0], data[:, 1], s=5, c=labels, cmap='Accent', alpha=0.5)  # 画出聚类后带标签的样本
plt.scatter(centers[:, 0], centers[:, 1], marker='*', s=120, c=[0, 1, 2], cmap='Accent')  # 画出聚类中心
plt.show()

2. 主成分分析PCA

PCA（Principal Component Analysis）是一种常见的数据分析方式，常用于高维数据的降维，可用于提取数据的主要特征分量。PCA 通常用于降低大型数据集的维数，使数据集中的指标数量变少，并且保留原数据集中指标的大部分信息。总而言之：减少数据指标数量，保留尽可能多的信息。

PCA 优点在于数据降维，便于提取数据的主要特征，使得数据更容易使用，减少计算开销，去除噪音等；PCA 适用于结构化数据，不仅能将数据压缩，也使得降维之后的数据特征相互独立。PCA 缺点在于不一定需要，有可能损失有用信息，只针对训练集保留主要信息，可能造成过拟合。

PCA 的步骤主要分为五步：标准化连续初始变量的范围、计算协方差矩阵以识别相关性、计算协方差矩阵的特征向量和特征值以识别主成分、创建特征向量来决定保留那些主成分、沿主成分轴重铸数据。

我们随机生成具有 y = 0.5x 的二维数据，并在 Y 轴方向上添加少量随机噪音。通过使用 PCA 法拟合数据，其中最重要的是分量和解释方差，对于这些数字的概念，让我们将其可视化为输入数据上的向量，使用分量来定义向量的方向，使用解释方差来定义向量的平方长度。

这些向量代表数据的主轴，向量的长度表明该轴在描述数据分布方面的重要性，更准确的说，它是投影时数据方差的度量到那个轴。每个数据点在主轴上的投影是数据的主成分。这种从数据轴到主轴的变换被称为 Affine transformation，基本上由平移，旋转和均匀缩放组成。

代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
def f(x):
    rs = np.random.RandomState(1)
    return x / 2 + rs.randn(*x.shape) / 4

rs = np.random.RandomState(1)
data = rs.randn(200, 2)
data[:, 1] = f(data[:, 0])
# plt.scatter(data[:, 0], data[:, 1], s=5, color='c')
# plt.axis('equal')
# plt.show()

pca = PCA(n_components=2)
pca.fit(data)
data_pca = pca.transform(data)
print(pca.components_)  # PCA分量
print(pca.explained_variance_)  # PCA解释方差

def draw_vector(v0, v1, ax=None):
    print(v0, v1)
    ax = ax or plt.gca()
    arrowprops = dict(arrowstyle='->', linewidth=2, shrinkA=0, shrinkB=0)
    ax.annotate('', v1, v0, arrowprops=arrowprops)

# 将分量和解释方差可视化为输入数据上的向量，使用分量来定义向量的方向，使用解释方差来定义向量的平方长度
plt.scatter(data[:, 0], data[:, 1], s=5, color='c', alpha=0.2)
for length, vector in zip(pca.explained_variance_, pca.components_):
    v = vector * 3 * np.sqrt(length)
    draw_vector(pca.mean_, pca.mean_ + v)
plt.axis('equal')
plt.show()

3. 高维数据PCA后聚类

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
# 产生高维数据
hd_centers = [
    [2, 1, 4, 2, 3, 2, 6, 1, 5, 3],
    [5, 7, 8, 5, 5, 9, 3, 6, 6, 5],
    [3, 8, 8, 9, 7, 1, 8, 8, 5, 4]
]
hd_center_num = len(hd_centers)
hd_data_num = 500

hd_data = np.zeros([hd_data_num * hd_center_num, 10])
for i, center in enumerate(hd_centers):
    hd_data[i * hd_data_num:(i + 1) * hd_data_num] = get_random_data(center, hd_data_num, data_dim=10)

print(hd_data.shape)  # (1500, 10)


# ---------- PCA主成分分析降维 ----------
pca = PCA(n_components=2)  # 降成二维
pca.fit(hd_data)
pca_data = pca.transform(hd_data)

print(pca_data.shape)  # (1500, 2)
plt.scatter(pca_data[:, 0], pca_data[:, 1], s=5, color='c')
plt.show()


# ---------- KMeans聚类 ----------
kms = KMeans(n_clusters=3, init='k-means++')
kms.fit(pca_data)  # 模型拟合
centers = kms.cluster_centers_  # 计算聚类中心
labels = kms.labels_  # 聚类后每个样本的类别标签

plt.figure(figsize=(8, 6))
plt.scatter(pca_data[:, 0], pca_data[:, 1], s=5, c=labels, cmap='Accent', alpha=0.5)  # 画出聚类后带标签的样本
plt.scatter(centers[:, 0], centers[:, 1], marker='*', s=120, c=[0, 1, 2], cmap='Accent')  # 画出聚类中心
plt.show()

4. 高维数据聚类并计算与中心的相似度

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
# 产生高维数据
hd_centers = [
    [2, 1, 4, 2, 3, 2, 6, 1, 5, 3],
    [5, 7, 8, 5, 5, 9, 3, 6, 6, 5],
    [3, 8, 8, 9, 7, 1, 8, 8, 5, 4]
]
hd_center_num = len(hd_centers)
hd_data_num = 500

hd_data = np.zeros([hd_data_num * hd_center_num, 10])
for i, center in enumerate(hd_centers):
    hd_data[i * hd_data_num:(i + 1) * hd_data_num] = get_random_data(center, hd_data_num, data_dim=10)

print(hd_data.shape)  # (1500, 10)


# 对原始高维数据聚类
kms = KMeans(n_clusters=3, init='k-means++')
kms.fit(hd_data)  # 模型拟合
centers = kms.cluster_centers_  # 计算聚类中心
labels = kms.labels_  # 聚类后每个样本的类别标签

print(centers.shape)  # (3, 10)，3个聚类中心的高维向量
print(labels[:10])  # [1 1 1 1 1 1 1 1 1 1]，第1类的标签
print(labels[500:510])  # [2 2 2 2 2 2 2 2 2 2]，第2类的标签
print(labels[1490:])  # [0 0 0 0 0 0 0 0 0 0]，第3类的标签


# 计算每个向量分别与3个聚类中心的余弦相似度
def CosineSimilarity(x, y):
    normalized_x = x / np.linalg.norm(x, axis=-1, keepdims=True)
    normalized_y = y / np.linalg.norm(y, axis=-1, keepdims=True)
    return np.sum(normalized_x * normalized_y, axis=-1)

data_cnt = hd_data_num * hd_center_num
true_cnt = 0

for i, data in enumerate(hd_data):
    score_0 = CosineSimilarity(data, centers[0])
    score_1 = CosineSimilarity(data, centers[1])
    score_2 = CosineSimilarity(data, centers[2])
    max_idx = np.argmax([score_0, score_1, score_2])
    print(f'ID: {i}, Score0: {score_0:.2f}, Score1: {score_1:.2f}, Score2: {score_2:.2f}, Pred_label: {max_idx}')
    if i < 500:
        if max_idx == 1:
            true_cnt = true_cnt + 1
    elif i < 1000:
        if max_idx == 2:
            true_cnt = true_cnt + 1
    else:
        if max_idx == 0:
            true_cnt = true_cnt + 1
# ID: 0, Score0: 0.85, Score1: 0.98, Score2: 0.85, Pred_label: 1
# ID: 1, Score0: 0.84, Score1: 0.92, Score2: 0.79, Pred_label: 1
# ID: 2, Score0: 0.87, Score1: 0.93, Score2: 0.90, Pred_label: 1
# ...
# ID: 500, Score0: 0.85, Score1: 0.78, Score2: 0.99, Pred_label: 2
# ID: 501, Score0: 0.85, Score1: 0.75, Score2: 0.98, Pred_label: 2
# ID: 502, Score0: 0.83, Score1: 0.71, Score2: 0.98, Pred_label: 2
# ...
# ID: 1000, Score0: 0.99, Score1: 0.81, Score2: 0.87, Pred_label: 0
# ID: 1001, Score0: 0.98, Score1: 0.77, Score2: 0.83, Pred_label: 0
# ID: 1002, Score0: 0.99, Score1: 0.76, Score2: 0.87, Pred_label: 0
# ...

print(f'All: {data_cnt}, True: {true_cnt}, Acc: {true_cnt / data_cnt:.2f}')  # All: 1500, True: 1499, Acc: 1.00