小小程序员
2021-05-13 15:56:57
机器学习之聚类算法
一、引言
聚类算法是无监督学习，只需要数据，而不需要标记结果，通过学习训练，用于发现共同的群体。本文将介绍几种常见的聚类算法，包括K-means、层次聚类和GMM高斯混合模型等。
1、聚类算法是什么
给定n个训练样本（未标记），如：{X1,X2,...,Xn}，同时给定聚类的个数K。
目标：把比较“接近”的样本放到一个簇类（cluster）里，总共得到K个簇类（cluster）。
2、如何衡量样本距离
下面是几个最常用的距离计算公式：
（1）欧氏距离：
（2）曼哈顿距离
（3）核函数映射后距离
3、如何评价聚类好坏
类间距高，类内距低，通俗地说就是“抱团紧不紧，异族远不远”。
4、应用场景
（1）图片检索：图片内容相似度
（2）图片分割：图片像素相似度
（3）网页聚类：文本内容相似度
（4）社交网络聚类：（被）关注人群、兴趣
（5）电商用户聚类：点击/加车/购买商品、行为序列等
二、K-means算法
K-means算法是使用最普遍，最重要的聚类算法之一。
已知观测集
1、算法原理
迭代收敛定义：
（1）聚类中心不在有变化
（2）每个样本到对应聚类中心的距离之和不再有很大变化
损失函数：
假定
为K个聚类中心，用
表示
是否属于聚类K，则损失函数定义如下：
最小化损失函数的过程是一个NP问题。上面的迭代，是收敛到局部最低点的过程。
2、优缺点
K-means算法优点：
（1）计算复杂度低，为O(Nmq)，其中N是数据总量，m是类别（即k），q是迭代次数
（2）原理简单，实现容易，容易解释
（3）聚类效果不错
K-means算法缺点：
（1）对异常值（噪声）敏感，可以通过一些调整（如中心值不直接取均值，而是找均值最近的样本点代替）
（2）需要提前确定K值（提前确定多少类）
（3）分类结果依赖于分类中心的初始化（可以通过进行多次K-means取最优来解决）
（4）属于“硬聚类”，每个样本只能有一个类别。其他的一些聚类方法（GMM或者模糊K-means允许“软聚类”）
（5）对于团状的数据点集区分度好，对于带状（环绕）等非凸形状不太好。（用谱聚类或做特征映射）
3、改进算法
如：K-means++等
4、动手用python实现K-means算法
样本数据：
1.658985 4.285136 -3.453687 3.424321 4.838138 -1.151539 -5.379713 -3.362104 0.972564 2.924086 -3.567919 1.531611 0.450614 -3.302219 -3.487105 -1.724432 2.668759 1.594842 -3.156485 3.191137 3.165506 -3.999838 -2.786837 -3.099354 4.208187 2.984927 -2.123337 2.943366 0.704199 -0.479481 -0.392370 -3.963704 2.831667 1.574018 -0.790153 3.343144 2.943496 -3.357075 -3.195883 -2.283926 2.336445 2.875106 -1.786345 2.554248 2.190101 -1.906020 -3.403367 -2.778288 1.778124 3.880832 -1.688346 2.230267 2.592976 -2.054368 -4.007257 -3.207066 2.257734 3.387564 -2.679011 0.785119 0.939512 -4.023563 -3.674424 -2.261084 2.046259 2.735279 -3.189470 1.780269 4.372646 -0.822248 -2.579316 -3.497576 1.889034 5.190400 -0.798747 2.185588 2.836520 -2.658556 -3.837877 -3.253815 2.096701 3.886007 -2.709034 2.923887 3.367037 -3.184789 -2.121479 -4.232586 2.329546 3.179764 -3.284816 3.273099 3.091414 -3.815232 -3.762093 -2.432191 3.542056 2.778832 -1.736822 4.241041 2.127073 -2.983680 -4.323818 -3.938116 3.792121 5.135768 -4.786473 3.358547 2.624081 -3.260715 -4.009299 -2.978115 2.493525 1.963710 -2.513661 2.642162 1.864375 -3.176309 -3.171184 -3.572452 2.894220 2.489128 -2.562539 2.884438 3.491078 -3.947487 -2.565729 -2.012114 3.332948 3.983102 -1.616805 3.573188 2.280615 -2.559444 -2.651229 -3.103198 2.321395 3.154987 -1.685703 2.939697 3.031012 -3.620252 -4.599622 -2.185829 4.196223 1.126677 -2.133863 3.093686 4.668892 -2.562705 -2.793241 -2.149706 2.884105 3.043438 -2.967647 2.848696 4.479332 -1.764772 -4.905566 -2.911070 |
其数据分布如下，可以看到，数据大概可以分成四类。
下面是通过K-means算法聚类后的分布图：
完整代码：k-means.py
# coding: utf-8
import random
import numpy as np
import matplotlib.pyplot as plt
def show_fig():
data = load_data()
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(data[:, 0], data[:, 1])
plt.show()
def cal_distance(node1, node2):
"""
计算两个向量之间的欧式距离
:param node1:
:param node2:
:return:
"""
return np.sqrt(np.sum(np.square(node1 - node2)))
def load_data():
"""
加载数据
:return:
"""
data = np.loadtxt("dataSet.csv")
return data
def init_k_node(data, k):
"""
随机初始化k个数据返回
:param data:
:param k:
:return:
"""
data = list(data)
return random.sample(data, k)
def get_clusters(data, k_centroids):
"""
计算各个节点所属的簇类
:param data:
:param k_centroids:
:return:
"""
cluster_dict = dict()
k = len(k_centroids)
for node in data:
# 计算节点node与k个质心的距离，选取距离最近的，并将节点加入相应簇类中
cluster_idx = -1
min_dis = float("inf")
for idx in range(k):
centroid = k_centroids[idx]
distance = cal_distance(node, centroid)
if distance < min_dis:
min_dis = distance
cluster_idx = idx
if cluster_idx not in cluster_dict.keys():
cluster_dict[cluster_idx] = []
cluster_dict[cluster_idx].append(node) # 加入对应类别中
return cluster_dict
def get_centroids(cluster_dict):
"""
重新计算k个质心
:param cluster_dict:
:return:
"""
new_k_centroids = []
for cluster_idx in cluster_dict.keys():
new_centroid = np.mean(cluster_dict[cluster_idx], axis=0)
new_k_centroids.append(new_centroid)
return new_k_centroids
def get_variance(centroids, cluster_dict):
"""
计算各个簇集合的均方误差
将簇类中各个节点与质心的距离累加求和
:param centroids:
:param cluster_dict:
:return:
"""
sum = 0.0
for cluster_idx in cluster_dict.keys():
centroid = centroids[cluster_idx]
distance = 0.0
for node in cluster_dict[cluster_idx]:
distance += cal_distance(node, centroid)
sum += distance
return sum
def show_cluster(centroids, cluster_dict):
"""
展示聚类结果
:param centroids:
:param cluster_dict:
:return:
"""
color_mark = ['or', 'ob', 'og', 'ok', 'oy', 'ow']
centroid_mark = ['dr', 'db', 'dg', 'dk', 'dy', 'dw']
for key in cluster_dict.keys():
plt.plot(centroids[key][0], centroids[key][1], centroid_mark[key], markersize=12) # 质心点
for node in cluster_dict[key]:
plt.plot(node[0], node[1], color_mark[key])
plt.show()
def main():
data = load_data()
centroids = init_k_node(data, 4)
cluster_dict = get_clusters(data, centroids)
new_var = get_variance(centroids, cluster_dict)
old_var = 1
# 当两次聚类的误差小于某个值时，说明质心基本稳定
while abs(new_var - old_var) >= 0.00001:
centroids = get_centroids(cluster_dict)
cluster_dict = get_clusters(data, centroids)
old_var = new_var
new_var = get_variance(centroids, cluster_dict)
show_cluster(centroids, cluster_dict)
if __name__ == '__main__':
show_fig()
main()
三、层次聚类
K-means一个最大的限制是，需要实现知道K值，即知道多少个分类。那么有没有不需要确定K值就可以分类的聚类算法呢？层次聚类就是这种算法。
层次聚类分为凝聚式层次聚类和分裂式层次聚类。
凝聚式层次聚类，就是在初始阶段将每一个点都视为一个簇，之后每一次合并两个最接近的簇。
分裂式层次聚类，就是在初始阶段将所有的点视为一个簇，之后每次分裂出一个簇，直到最后剩下单个点的簇为止。
1、簇之间的距离定义
下面讲解凝聚式聚类。凝聚式聚类簇之间的距离定义有如下三种。
cluster R和cluster S之间的距离计算有以下几种：
（1）最小连接距离法（Single linkage）
取两个簇中距离最近的两个样本的距离作为这两个簇的距离。
（2）最大连接距离发（Complete linkage）
取两个簇中距离最远的两个点的距离作为这两个簇的距离。
（3）平均连接距离法（Average linkage）
把两个簇中的点的两两的距离全部放在一起求均值，取得的结果也相对合适一点。
2、层次聚类算法原理
（1）初始化，将每个样本都视为一个簇
（2）计算各个簇之间的距离
（3）寻找最近的两个簇，并将他们合并为一类
（4）重复步骤（2）和步骤（3），知道所有样本归为一类。
3、层次聚类算法过程举例
如下图是5个点的样本分布图：
初始化，将每个样本视为一个簇，其欧式距离原始举证如下：
| P1 | P2 | P3 | P4 | P5 | |
| P1 | 0 | 0.81 | 1.32 | 1.94 | 1.82 |
| P2 | 0.81 | 0 | 1.56 | 2.16 | 1.77 |
| P3 | 1.32 | 1.56 | 0 | 0.63 | 0.71 |
| P4 | 1.94 | 2.16 | 0.63 | 0 | 0.71 |
| P5 | 1.82 | 1.77 | 0.71 | 0.71 | 0 |
按照算法的流程，我们取距离最近的两个簇P3和P4并合并为{P3,P4}，根据最小连接距离法更新距离矩阵：
| P1 | P2 | {P3,P4} | P5 | |
| P1 | 0 | 0.81 | 1.32 | 1.82 |
| P2 | 0.81 | 0 | 1.56 | 1.77 |
| {P3,P4} | 1.32 | 1.56 | 0 | 0.71 |
| P5 | 1.82 | 1.77 | 0.71 | 0 |
接着找出距离最近的两个簇{P3,P4}和P5，合并为{P3,P4,P5}，根据最小连接距离法更新距离矩阵：
| P1 | P2 | {P3,P4,P5} | |
| P1 | 0 | 0.81 | 1.32 |
| P2 | 0.81 | 0 | 1.56 |
| {P3,P4,P5} | 1.32 | 1.56 | 0 |
接着找出距离最近的两个簇P1和P2，合并为{P1,P2}，根据最小连接距离法更新距离矩阵：
| {P1,P2} | {P3,P4,P5} | |
| {P1,P2} | 0 | 1.32 |
| {P3,P4,P5} | 1.32 | 0 |
最终合并剩下的两个簇，得到如下图的聚类结果：
4、用python实现层次聚类算法
#-*- coding:utf-8 -*-
import math
import pylab as pl
def dist(node1, node2):
"""
计算欧几里得距离,node1,node2分别为两个元组
:param node1:
:param node2:
:return:
"""
return math.sqrt(math.pow(node1[0]-node2[0], 2)+math.pow(node1[1]-node2[1], 2))
def dist_min(cluster_x, cluster_y):
"""
Single linkage
又叫做 nearest-neighbor ，就是取两个类中距离最近的两个样本的距离作为这两个集合的距离。
:param cluster_x:
:param cluster_y:
:return:
"""
return min(dist(node1, node2) for node1 in cluster_x for node2 in cluster_y)
def dist_max(cluster_x, cluster_y):
"""
Complete linkage
这个则完全是 Single linkage 的反面极端，取两个集合中距离最远的两个点的距离作为两个集合的距离。
:param cluster_x:
:param cluster_y:
:return:
"""
return max(dist(node1, node2) for node1 in cluster_x for node2 in cluster_y)
def dist_avg(cluster_x, cluster_y):
"""
Average linkage
这种方法就是把两个集合中的点两两的距离全部放在一起求均值，相对也能得到合适一点的结果。
:param cluster_x:
:param cluster_y:
:return:
"""
return sum(dist(node1, node2) for node1 in cluster_x for node2 in cluster_y)/(len(cluster_x)*len(cluster_y))
def find_min(distance_matrix):
"""
找出距离最近的两个簇下标
:param distance_matrix:
:return:
"""
min = 1000
x = 0; y = 0
for i in range(len(distance_matrix)):
for j in range(len(distance_matrix[i])):
if i != j and distance_matrix[i][j] < min:
min = distance_matrix[i][j];x = i; y = j
return (x, y, min)
def AGNES(dataset, distance_method, k):
"""
聚类算法模型
:param dataset: 数据集
:param distance_method: 计算簇类聚类的方法
:param k: 目标簇类数目
:return:
"""
print len(dataset)
# 初始化簇类集合和距离矩阵
cluster_set = []
distance_matrix = []
for node in dataset:
cluster_set.append([node])
print 'original cluster set:'
print cluster_set
for cluster_x in cluster_set:
distance_list = []
for cluster_y in cluster_set:
distance_list.append(distance_method(cluster_x, cluster_y))
distance_matrix.append(distance_list)
print 'original distance matrix:'
print len(distance_matrix), len(distance_matrix[0]), distance_matrix
q = len(dataset)
# 合并更新
while q > k:
id_x, id_y, min_distance = find_min(distance_matrix)
cluster_set[id_x].extend(cluster_set[id_y])
cluster_set.remove(cluster_set[id_y])
distance_matrix = []
for cluster_x in cluster_set:
distance_list = []
for cluster_y in cluster_set:
distance_list.append(distance_method(cluster_x, cluster_y))
distance_matrix.append(distance_list)
q -= 1
return cluster_set
def draw(cluster_set):
"""
画图
:param cluster_set:
:return:
"""
color_list = ['r', 'y', 'g', 'b', 'c', 'k', 'm']
for cluster_idx, cluster in enumerate(cluster_set):
coo_x = [] # x坐标列表
coo_y = [] # y坐标列表
for node in cluster:
coo_x.append(node[0])
coo_y.append(node[1])
pl.scatter(coo_x, coo_y, marker='x', color=color_list[cluster_idx % len(color_list)], label=cluster_idx)
pl.legend(loc='upper right')
pl.show()
# 数据集：每三个是一组分别是西瓜的编号，密度，含糖量
data = """
1,0.697,0.46,2,0.774,0.376,3,0.634,0.264,4,0.608,0.318,5,0.556,0.215,
6,0.403,0.237,7,0.481,0.149,8,0.437,0.211,9,0.666,0.091,10,0.243,0.267,
11,0.245,0.057,12,0.343,0.099,13,0.639,0.161,14,0.657,0.198,15,0.36,0.37,
16,0.593,0.042,17,0.719,0.103,18,0.359,0.188,19,0.339,0.241,20,0.282,0.257,
21,0.748,0.232,22,0.714,0.346,23,0.483,0.312,24,0.478,0.437,25,0.525,0.369,
26,0.751,0.489,27,0.532,0.472,28,0.473,0.376,29,0.725,0.445,30,0.446,0.459"""
# 数据处理 dataset是30个样本（密度，含糖量）的列表
a = data.split(',')
dataset = [(float(a[i]), float(a[i+1])) for i in range(1, len(a)-1, 3)]
cluster_set = AGNES(dataset, dist_avg, 3)
print 'final cluster set:'
print cluster_set
draw(cluster_set)最终效果：
四、GMM高斯混合模型和EM算法
EM算法是一种迭代算法，1977年由Dempster等人总结提出，用于含有隐变量的概率模型参数的极大似然估计，或极大后验概率估计。
EM算法的每次迭代由两部组成：E步，求期望；M步，求极大。所以这一算法称为期望极大算法，简称EM算法。
4.1 高斯混合模型
定义：高斯混合模型是指具有如下形式的概率分布模型：
其中，
是系数，大于等于0，
；
是高斯分布密度，
，
称为第k个分模型。
4.2 高斯混合模型参数估计的EM算法
输入：假设观测数据
，由高斯混合模型生成，即：
其中
输出：利用EM算法估计高斯混合模型的参数θ。
（1）取参数的初始值开始迭代
（2）E步：依据当前模型参数，计算分模型k对观测数据
的响应度
（3）M步：计算新一轮的模型参数
（4）重复（2）步和（3）步，直到收敛。
4.3 python实现
生成样本数据：
import numpy as np
import matplotlib.pyplot as plt
cov1 = np.mat("0.3 0;0 0.1")
cov2 = np.mat("0.2 0;0 0.3")
mu1 = np.array([0, 1])
mu2 = np.array([2, 1])
sample = np.zeros((100, 2))
sample[:30, :] = np.random.multivariate_normal(mean=mu1, cov=cov1, size=30)
sample[30:, :] = np.random.multivariate_normal(mean=mu2, cov=cov2, size=70)
np.savetxt("sample.data", sample)
plt.plot(sample[:30, 0], sample[:30, 1], "bo")
plt.plot(sample[30:, 0], sample[30:, 1], "rs")
plt.show()
原始数据分布图
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
######################################################
# 第 k 个模型的高斯分布密度函数
# 每 i 行表示第 i 个样本在各模型中的出现概率
# 返回一维列表
######################################################
def phi(Y, mu_k, cov_k):
norm = multivariate_normal(mean=mu_k, cov=cov_k)
return norm.pdf(Y)
######################################################
# E 步：计算每个模型对样本的响应度
# Y 为样本矩阵，每个样本一行，只有一个特征时为列向量
# mu 为均值多维数组，每行表示一个样本各个特征的均值
# cov 为协方差矩阵的数组
# alpha 为模型响应度数组
######################################################
def getExpectation(Y, mu, cov, alpha):
# 样本数
N = Y.shape[0]
# 模型数
K = alpha.shape[0]
# 为避免使用单个高斯模型或样本，导致返回结果的类型不一致
# 因此要求样本数和模型个数必须大于1
assert N > 1, "There must be more than one sample!"
assert K > 1, "There must be more than one gaussian model!"
# 响应度矩阵，行对应样本，列对应响应度
gamma = np.mat(np.zeros((N, K)))
# 计算各模型中所有样本出现的概率，行对应样本，列对应模型
prob = np.zeros((N, K))
for k in range(K):
prob[:, k] = phi(Y, mu[k], cov[k])
prob = np.mat(prob)
# 计算每个模型对每个样本的响应度
for k in range(K):
gamma[:, k] = alpha[k] * prob[:, k]
for i in range(N):
gamma[i, :] /= np.sum(gamma[i, :])
return gamma
######################################################
# M 步：迭代模型参数
# Y 为样本矩阵，gamma 为响应度矩阵
######################################################
def maximize(Y, gamma):
# 样本数和特征数
N, D = Y.shape
# 模型数
K = gamma.shape[1]
#初始化参数值
mu = np.zeros((K, D))
cov = []
alpha = np.zeros(K)
# 更新每个模型的参数
for k in range(K):
# 第 k 个模型对所有样本的响应度之和
Nk = np.sum(gamma[:, k])
# 更新 mu
# 对每个特征求均值
for d in range(D):
mu[k, d] = np.sum(np.multiply(gamma[:, k], Y[:, d])) / Nk
# 更新 cov
cov_k = np.mat(np.zeros((D, D)))
for i in range(N):
cov_k += gamma[i, k] * (Y[i] - mu[k]).T * (Y[i] - mu[k]) / Nk
cov.append(cov_k)
# 更新 alpha
alpha[k] = Nk / N
cov = np.array(cov)
return mu, cov, alpha
######################################################
# 数据预处理
# 将所有数据都缩放到 0 和 1 之间
######################################################
def scale_data(Y):
# 对每一维特征分别进行缩放
for i in range(Y.shape[1]):
max_ = Y[:, i].max()
min_ = Y[:, i].min()
Y[:, i] = (Y[:, i] - min_) / (max_ - min_)
print("Data scaled.")
return Y
######################################################
# 初始化模型参数
# shape 是表示样本规模的二元组，(样本数, 特征数)
# K 表示模型个数
######################################################
def init_params(shape, K):
N, D = shape
mu = np.random.rand(K, D) # 生成K个，D维的随机样本数据，并且随机样本位于[0,1)中
cov = np.array([np.eye(D)] * K) # 生成包含K个DxD对角矩阵的矩阵
alpha = np.array([1.0 / K] * K) # 生成K个值为1/K的矩阵
# print "Parameters initialized."
# print "mu:", mu
# print "cov:", cov
# print "alpha:", alpha
return mu, cov, alpha
######################################################
# 高斯混合模型 EM 算法
# 给定样本矩阵 Y，计算模型参数
# K 为模型个数
# times 为迭代次数
######################################################
def GMM_EM(Y, K, times):
Y = scale_data(Y)
mu, cov, alpha = init_params(Y.shape, K)
for i in range(times):
gamma = getExpectation(Y, mu, cov, alpha)
mu, cov, alpha = maximize(Y, gamma)
print("{sep} Result {sep}".format(sep="-" * 20))
print "mu:", mu
print "cov:", cov
print "alpha:", alpha
return mu, cov, alpha
# 载入数据
Y = np.loadtxt("gmm.data")
matY = np.matrix(Y, copy=True)
# 模型个数，即聚类的类别个数
K = 2
# 计算 GMM 模型参数
mu, cov, alpha = GMM_EM(matY, K, 100)
# 根据 GMM 模型，对样本数据进行聚类，一个模型对应一个类别
N = Y.shape[0]
# 求当前模型参数下，各模型对样本的响应度矩阵
gamma = getExpectation(matY, mu, cov, alpha)
# 对每个样本，求响应度最大的模型下标，作为其类别标识
category = gamma.argmax(axis=1).flatten().tolist()[0]
# 将每个样本放入对应类别的列表中
class1 = np.array([Y[i] for i in range(N) if category[i] == 0])
class2 = np.array([Y[i] for i in range(N) if category[i] == 1])
# 绘制聚类结果
plt.plot(class1[:, 0], class1[:, 1], 'rs', label="class1")
plt.plot(class2[:, 0], class2[:, 1], 'bo', label="class2")
plt.legend(loc="best")
plt.title("GMM Clustering By EM Algorithm")
plt.show()
聚类结果分布图
对比原数据的分布图和聚类的分布图，可以看到效果还是非常好的。
条内容