一、引言

聚类算法是无监督学习,只需要数据,而不需要标记结果,通过学习训练,用于发现共同的群体。本文将介绍几种常见的聚类算法,包括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算法是使用最普遍,最重要的聚类算法之一。

已知观测集,其中每个观测都是一个

D

-维实向量,k-平均聚类要把这个观测划分到k个集合中(k≤n),使得组内平方和(WCSS within-cluster sum of squares)最小。换句话说,它的目标是找到使得下式满足的聚类

{\underset {{\mathbf {S}}}{\operatorname {arg\,min}}}\sum _{{i=1}}^{{k}}\sum _{{{\mathbf x}\in S_{i}}}\left\|{\mathbf x}-{\boldsymbol \mu }_{i}\right\|^{2}

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. 层次聚类计算准则

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()

聚类结果分布图

对比原数据的分布图和聚类的分布图,可以看到效果还是非常好的。

 

 

打赏

发表评论

电子邮件地址不会被公开。