Clustering杂谈

Posted by SkyHigh on October 22, 2016

Clustering杂谈

聚类算法有很多,主要的目标就是将一堆数据中相似的数据归到一起,而不去管这个数据的具体标签。归类的方法直接从数据点中获取,在机器学习上这种方法属于非监督学习(Unsupervised Learning)。
常用的聚类算法,比较经典的有k-means、k-medoids、GMM等。比较现代和流行的是Spectral Clustering(谱聚类)和Hierarchical Clustering(层级聚类)。

本节用的例子数据地址:

pickle数据

直接复制粘贴即可。

所有代码地址:

Clustering代码

1.k-means

k均值聚类的思路比较清晰,就是给每个数据分配一个具体的标签。使用k均值聚类需要事先确定聚类数。
k-means的算法步骤如下:

  1. 选定$K$个中心$\mu_k$的初值。这个过程通常是针对具体的问题有一些启发式的选取方法,或者大多数情况下采用随机选取的办法。因为前面说过k-means并不能保证全局最优,而是否能收敛到全局最优解其实和初值的选取有很大的关系,所以有时候我们会多次选取初值跑k-means,并取其中最好的一次结果。
  2. 将每个数据点归类到离它最近的那个中心点所代表的cluster中。 用公式$\mu_k = \frac{1}{N_k}\sum_{j\in\text{cluster}_k}x_j $计算出每个cluster的新的中心点。
  3. 重复第二步,一直到迭代了最大的步数或者前后的总误差$J$的值相差小于一个阈值为止。其中总误差$J$的公式如下:
    • $\displaystyle J = \sum_{n=1}^N\sum_{k=1}^K r_{nk} |x_n-\mu_k|^2$

k-means的迭代过程其实是简化的EM算法思想。

代码

# -*- encoding:utf-8 -*-
import numpy as np
from matplotlib import pyplot as plt
import pickle
from copy import deepcopy
import random

class Kmeans():
    def __init__(self, k=3, filepath='../data/cluster.pkl'):
        self.samples = self.loadData(filepath)
        # K个聚类点component
        self.K = k
        # 数据预处理
        self.random_points = np.zeros([len(self.samples)*len(self.samples[0][0]), 2])
        i = 0
        for smp in self.samples:
            smp_len = len(smp[0])
            self.random_points[i:i+smp_len,0] = smp[0]
            self.random_points[i:i+smp_len,1] = smp[1]
            i += smp_len
        # M条数据,每条数据N个特征
        self.M, self.N = self.random_points.shape
        self.centroids = self.initial(self.K)
        self.dist = np.zeros((self.M, self.K))

    # 加载数据
    def loadData(self, filepath):
        with open(filepath) as fr:
            return pickle.load(fr)

    # k均值
    def kmeans(self, K, all=None):
    	# 初始化中心点
        centroids = self.initial(K, all)
        if not all:
            all = range(0, self.M)
        dist = deepcopy(self.dist)
        cost = np.inf
        # 循环终止门限
        threshold = 1e-3
        while True:
            cur_cost = 0.
            # 计算每个点到当前各个聚类点的距离
            for k in range(K):
                dist[all,k] = np.sqrt(np.sum(np.power(self.random_points[all,:]-centroids[k], 2), 1))
            # 更新中心点
            for k in range(K):
                all_ind = np.argmin(dist, 1)
                k_ind = (all_ind == k).nonzero()[0]
                k_ind = list((set(all) & set(k_ind)))
                k_points = self.random_points[k_ind,:]
                centroids[k] = np.average(k_points, 0)
                cur_cost += np.sum(dist[k_ind, k])

            print cur_cost
            # 如果收敛,则跳出循环
            if abs(cur_cost-cost)<threshold:
                break

            cost = cur_cost

        # self.centroids = centroids
        # self.dist = dist
        return dist, centroids

    def initial(self, K, all=None):
        if not all:
            return self.random_points[np.random.choice(range(0,self.M), K, replace=False)]
        else:
            return self.random_points[np.random.choice(all, K, replace=False)]

    def draw(self, dist):
        plt.figure(1)
        all_ind = np.argmin(dist, 1)
        color = ['r', 'g', 'b']
        for k in range(self.K):
            x, y = self.random_points[(all_ind==k).nonzero()[0],0],self.random_points[(all_ind==k).nonzero()[0],1]
            plt.scatter(x, y , c=color[k])
        plt.show()


if __name__ == '__main__':
    km = Kmeans()
    dist, centroids = km.kmeans(3)
    km.draw(dist)

结果

原数据

聚类结果

二分k-means

关于k-means的变种,还有二分k-means。二分k-means的算法步骤如下:

  1. 将所有数据点看成一个簇(初始化的簇位置为所有点的均值)。
  2. 对每一个簇,计算该簇的总误差,并对该簇进行k-均值聚类(k=2),计算将该簇一分为二后的总误差和。误差可以是簇中所有点到簇的聚类点的平均距离
  3. 选择使得误差最小的那个簇进行划分操作。
  4. 循环2到3,知道簇的数目为k为止。

二分k-means代码

# -*- encoding:utf-8 -*-
import numpy as np
from matplotlib import pyplot as plt
import pickle
from copy import deepcopy
import random

class Kmeans():
    def __init__(self, k=3, filepath='../data/cluster.pkl'):
        self.samples = self.loadData(filepath)
        # K个聚类点component
        self.K = k
        # 数据预处理
        self.random_points = np.zeros([len(self.samples)*len(self.samples[0][0]), 2])
        i = 0
        for smp in self.samples:
            smp_len = len(smp[0])
            self.random_points[i:i+smp_len,0] = smp[0]
            self.random_points[i:i+smp_len,1] = smp[1]
            i += smp_len
        # M条数据,每条数据N个特征
        self.M, self.N = self.random_points.shape

    # 加载数据
    def loadData(self, filepath):
        with open(filepath) as fr:
            return pickle.load(fr)

    # k均值
    def kmeans(self, K, all=None):
        centroids = self.initial(K, all)
        if all is None:
            all = range(0, self.M)
        all = np.array(all)
        dist = np.zeros([len(all), K])
        cost = np.inf
        threshold = 1e-3
        while True:
            cur_cost = 0.
            split_parts = list()
            split_cost = list()
            split_dist_ind = list()
            for k in range(K):
                dist[:,k] = np.sqrt(np.sum(np.power(self.random_points[all,:]-centroids[k], 2), 1))
            for k in range(K):
                all_ind = np.argmin(dist, 1)
                dist_ind = (all_ind==k).nonzero()[0]
                k_ind = all[dist_ind]
                split_parts.append(k_ind)
                k_points = self.random_points[k_ind,:]
                centroids[k] = np.average(k_points, 0)
                cur_cost += np.sum(dist[dist_ind, k])
                split_cost.append(np.sum(dist[dist_ind, k]))
                split_dist_ind.append(dist_ind)

            print cur_cost
            if abs(cur_cost-cost)<threshold:
                break

            cost = cur_cost
            del split_parts
            del split_cost

        # self.centroids = centroids
        # self.dist = dist
        return dist, centroids, split_parts, cur_cost, split_cost


    def bikmeans(self, K):
        # 簇的数量
        cur_num = 1
        # 各个簇的聚类点
        centroids = [list(np.average(self.random_points, 0))]
        # 各个簇的cost
        cost = [np.sum(np.sqrt(np.sum(np.power(self.random_points-centroids[0], 2), 1)))]
        # 各个簇的数据样本ID
        all = [range(0, self.M)]
        # 对于K个聚类点,进行K-1次划分
        for _ in range(K-1):
            min_cost = np.inf
            # 对每个簇进行遍历更新
            for i in range(cur_num):
                # 当前簇的聚类点(2个),划分的数据,划分前的cost,划分后的cost
                _, cur_centroids, split_parts, cur_cost, split_cost = self.kmeans(2, all[i])
                if i == 0 or cost[i]-cur_cost > cost[min_ind]-min_cost:
                    min_cost = cur_cost
                    min_ind = i
                    min_centroids = cur_centroids
                    min_split_cost = split_cost
                    min_split_parts = split_parts
            centroids[min_ind] = min_centroids[0]
            centroids.append(min_centroids[1])
            all[min_ind] = min_split_parts[0]
            all.append(min_split_parts[1])
            cost[min_ind] = min_split_cost[0]
            cost.append(min_split_cost[1])

            cur_num += 1

        return all
    
    
    def initial(self, K, all=None):
        if all is None:
            return self.random_points[np.random.choice(range(0,self.M), K, replace=False)]
        else:
            return self.random_points[np.random.choice(all, K, replace=False)]

    def drawBiKmeans(self, points):
        plt.figure(1)
        color = ['r', 'g', 'b']
        for k in range(self.K):
            x, y = self.random_points[points[k],0],self.random_points[points[k],1]
            plt.scatter(x, y , c=color[k])
        plt.show()


if __name__ == '__main__':
    km = Kmeans()
    points = km.bikmeans(3)
    km.drawBiKmeans(points)

二分kmeans结果

2.k-medoids

k均值的簇坐标为所有簇点的均值(means),而k-medoids则是选取簇中存在的点,该点是距离同个簇的所有其他点距离和最小的点,即寻找中位数(medium)。
k-medoids对于聚类点的选取要求没有k-means高,因此k-means能做的k-medoids也能做。k-medoids主要针对的是非数字化特征的数据,每个样本点之间只有一个差异度值。那么我们就可以构造一个距离矩阵(N*K),相应的误差公式如下:

$\mathcal{V}(x_n,\mu_k)$表示数据点到中心点的距离。而$r_{nk}$表示该数据点n是否在中心点k的簇里,是为1,否则为0。

可以看到k-means只要求得簇点均值即可,复杂度为O(N)。而k-medoids则要计算簇中每个点到其他点的距离,复杂度为O(N^2)。但是由于k-medoids找的是已有的点,它不会像k-means一样产生离群(outlier)的聚类点,也不会选择数据中的离群点为聚类点,因此k-medoids比k-means要更具健壮(robust)。

代码

# -*- encoding:utf-8 -*-
import numpy as np
from matplotlib import pyplot as plt
import pickle
from copy import deepcopy
import random

class Kmeans():
    def __init__(self, k=3, filepath='../data/cluster.pkl'):
        self.samples = self.loadData(filepath)
        # K个聚类点component
        self.K = k
        # 数据预处理
        self.random_points = np.zeros([len(self.samples)*len(self.samples[0][0]), 2])
        i = 0
        for smp in self.samples:
            smp_len = len(smp[0])
            self.random_points[i:i+smp_len,0] = smp[0]
            self.random_points[i:i+smp_len,1] = smp[1]
            i += smp_len
        # M条数据,每条数据N个特征
        self.M, self.N = self.random_points.shape

    # 加载数据
    def loadData(self, filepath):
        with open(filepath) as fr:
            return pickle.load(fr)

    def kmedoids(self, K):
        # 初始化聚类点
        centroids = self.initial(K)
        dist = np.zeros([self.M, K])
        cost = np.inf
        threshold = 1e-3
        while True:
            cur_cost = 0.
            # 计算每个点到当前各个聚类点的距离
            for k in range(K):
                dist[:,k] = np.sqrt(np.sum(np.power(self.random_points-centroids[k], 2), 1))
            # 更新聚类点
            for k in range(K):
                all_ind = np.argmin(dist, 1)
                dist_ind = (all_ind==k).nonzero()[0]
                k_points = self.random_points[dist_ind,:]
                dist_sum = np.zeros([len(k_points)])
                # 这里和kmeans不同,计算每个点到其他点的距离
                for i, point in enumerate(k_points):
                    dist_sum[i] = np.sum(np.sqrt(np.sum(np.power((self.random_points[dist_ind,:]-point), 2),1)))
                # 选出距离和最小的点作为聚类点
                best = dist_sum.argmin()
                centroids[k] = k_points[best]
                cur_cost += np.sum(dist_sum[best])

            print cur_cost
            # 如果收敛,则跳出循环
            if abs(cur_cost-cost)<threshold:
                break

            cost = cur_cost

        # self.centroids = centroids
        # self.dist = dist
        return dist, centroids


    def initial(self, K, all=None):
        if all is None:
            return self.random_points[np.random.choice(range(0,self.M), K, replace=False)]
        else:
            return self.random_points[np.random.choice(all, K, replace=False)]

    def draw(self, dist):
        plt.figure(1)
        all_ind = np.argmin(dist, 1)
        color = ['r', 'g', 'b']
        for k in range(self.K):
            x, y = self.random_points[(all_ind==k).nonzero()[0],0],self.random_points[(all_ind==k).nonzero()[0],1]
            plt.scatter(x, y , c=color[k])
        plt.show()

if __name__ == '__main__':
    km = Kmeans()
    dist, _ = km.kmedoids(3)
    km.draw(dist)

结果

原数据

聚类结果

3.Gaussian Mixture Model(GMM)

k-means是给每个数据点分配到具体的簇,而GMM的思路则是更贝叶斯。它不给出每个数据点对应的簇,而是给出每个数据点由各个聚类点产生的概率,并且不同聚类点产生某个数据点的概率和为1。
GMM的目的是,假定k个簇中每个簇产生的数据点满足n维高斯分布,估计产生目前观测的数据点的概率最大的k个簇的参数。即我们对下面的公式求最大值:

考虑到各个数据点的分布独立,则有:

考虑单个数据点的概率比较小,许多很小的数字相乘起来在计算机里很容易造成浮点数下溢,因此我们使用对数处理上式:

下面通过EM算法迭代求解参数$\mu_k$、$\pi_k$和$\Sigma_k$:

  1. E步
    • 对于每个数据$x_i$来说,它由第k个聚类点生成的概率为:
    • 在计算$\gamma(i, k)$的时候我们假定$\mu_k$和$\Sigma_k$均已知,我们将取上一次迭代所得的值(或者初始值)
  2. M步
    • 其中$N_k = \sum_{i=1}^N \gamma(i, k)$,N为总数据点数。

重复迭代E、M两步,直到对数似然函数的值收敛为止。一般到M步停止收敛。

当$\Sigma$为奇异矩阵时,即它的行列式等于0,此时计算概率的时候会出现问题(概率会等于0),需要另外处理。比如可以在$\Sigma$的对角线上加一个较小值。具体可以看参考链接。

实现

# -*- encoding:utf-8 -*-
import numpy as np
from numpy import linalg
import sys
import time
import pickle
reload(sys)
sys.setdefaultencoding('utf8')


class GMM():
    def __init__(self, k=3, filepath='../data/cluster.pkl'):
        self.samples = self.loadData(filepath)
        # K个聚类点component
        self.K = k
        # 数据预处理
        self.random_points = np.zeros([len(self.samples)*len(self.samples[0][0]), 2])
        i = 0
        for smp in self.samples:
            smp_len = len(smp[0])
            self.random_points[i:i+smp_len,0] = smp[0]
            self.random_points[i:i+smp_len,1] = smp[1]
            i += smp_len
        # M条数据,每条数据N个特征
        self.M, self.N = self.random_points.shape

    # 用k均值的思路初始化所有参数
    def initial(self):
        centroids = np.random.choice(range(self.M), self.K)
        centroids = self.random_points[centroids]

        # 1*K,即各个聚类的概率
        ppi = np.zeros([self.K])
        # K*N,即各个聚类点坐标
        pmiu = centroids
        # K*N*N,即各个聚类点的方差
        psigma = np.zeros([self.K, self.N, self.N])

        # M*K,距离矩阵
        distmat = np.tile(np.sum(self.random_points**2, 1).reshape([-1,1]), (1, self.K)) + \
                    np.tile(np.sum(pmiu**2, 1).reshape([1,-1]), (self.M, 1)) - \
                    2*self.random_points.dot(np.transpose(pmiu))

        # 初始化参数
        distind = distmat.argmin(1)

        for k in range(self.K):
            # 得到每个类的索引
            cur_ind = (distind==k).nonzero()[0]
            # 得到每个类的数据比例
            ppi[k] = float(len(cur_ind))/len(distind)
            # 得到每个类的协方差矩阵,N*N
            psigma[k,:,:] = np.cov(self.random_points[cur_ind].transpose())

        return ppi, pmiu, psigma



    def gmm(self):
        ppi, pmiu, psigma = self.initial()

        # 门限
        threshold = 1e-10

        Loss = -np.inf

        iter = 1
        while True:
            # M*K
            Pc = self.getProb(ppi, pmiu, psigma)

            # M*K
            gama = Pc*ppi
            gama = gama / np.sum(gama, 1).reshape(-1,1)

            # K*1
            Nk = np.sum(gama, 0).flatten()


            # 先计算协方差矩阵
            for k in range(self.K):
                shift = self.random_points-pmiu[k,:]
                psigma[k,:,:] = (np.transpose(shift).dot(np.diag(gama[:,k]).dot(shift)))/Nk[k]

            # 计算pi
            ppi = Nk/self.M
            # 计算miu
            pmiu = (np.dot(np.transpose(gama), self.random_points))/Nk.reshape([-1,1])

            # 也不一定是损失,只要使用一个会收敛的数值进行收敛判断就行,比如各个点到各自聚类点的距离之和
            curLoss = -np.sum(np.log(Pc.dot(ppi.reshape(self.K, 1))))

            # 判定是否收敛
            if abs(curLoss-Loss)<threshold:
                break
            Loss = curLoss
            print "Cur Loss:",Loss
            iter += 1

        self.Pc = Pc
        self.miu = pmiu
        self.gama = gama
        self.pi = ppi

        return Pc

    def getProb(self, ppi, pmiu, psigma):
        # N(x|pMiu,pSigma) = 1/((2pi)^(N/2))*(1/(abs(sigma))^0.5)*exp(-1/2*(x-pmiu)'psigma^(-1)*(x-pmiu))
        Pc = np.zeros([self.M, self.K])
        pi = np.pi
        N = self.N

        # 根据n维高斯公式计算概率密度矩阵
        for k in range(self.K):
            curp = psigma[k,:,:]
            sigma_inverse = linalg.pinv(curp)
            sigma_det = linalg.det(sigma_inverse)
            if sigma_det < 0:
                sigma_det = 0.
            shift = self.random_points - pmiu[k,:]
            Pc[:,k] = (1./((2*pi)**(N/2)))*(np.sqrt(np.abs(sigma_det)))*np.exp(-0.5*np.sum(np.dot(shift,sigma_inverse)*shift, 1)).flatten()
        # M*K
        return Pc

    # 加载数据
    def loadData(self, filepath):
        with open(filepath) as fr:
            return pickle.load(fr)

    # 绘制图形
    def draw(self, dataxy):
        from matplotlib import pyplot as plt
        plt.figure(1)
        plt.scatter(dataxy[:,0],dataxy[:,1])
        plt.show()

    # 绘制聚类图形
    def drawCluster(self):
        from matplotlib import pyplot as plt
        plt.figure(1)
        all_ind = np.argmax(self.Pc, 1)
        color = ['r', 'g', 'b']
        for k in range(self.K):
            print (all_ind==k).nonzero()[0]
            x, y = self.random_points[(all_ind==k).nonzero()[0],0],self.random_points[(all_ind==k).nonzero()[0],1]
            plt.scatter(x, y , c=color[k])
        plt.show()




if __name__ == '__main__':
    gmm = GMM()
    gmm.gmm()
    gmm.draw(gmm.random_points)
    gmm.drawCluster()

结果

原数据

聚类结果

参考