<commit_before>--- 
layout: post 
title: Clustering杂谈
date: 2016-10-21 
categories: blog 
tags: [聚类, 机器学习] 
description: 聚类算法
--- 

# Clustering杂谈

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

### 二分k-means
 
关于k-means的变种，还有二分k-means。二分k-means的算法步骤如下：


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

## 2.k-medoids

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

$$\displaystyle\tilde{J} = \sum_{n=1}^N\sum_{k=1}^K r_{nk}\mathcal{V}(x_n,\mu_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）。

## 3.Gaussian Mixture Model（GMM）

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

$$\displaystyle
\begin{aligned}
p(x) & = \sum_{k=1}^K p(k)p(x|k) \\
& = \sum_{k=1}^K \pi_k \mathcal{N}(x|\mu_k, \Sigma_k)
\end{aligned}
$$

考虑到各个数据点的分布独立，则有：

$$p(x)=\prod_{i=1}^N p(x_i)=\displaystyle
\sum_{i=1}^N \sum_{k=1}^K \pi_k \mathcal{N}(x_i|\mu_k, \Sigma_k)
$$

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

$$\log p(x)=\displaystyle
\sum_{i=1}^N \log \left\{\sum_{k=1}^K \pi_k \mathcal{N}(x_i|\mu_k, \Sigma_k)\right\}
$$

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

1. E步
	* 对于每个数据$x_i$来说，它由第k个聚类点生成的概率为：
	$$\displaystyle
	\gamma(i, k) = \frac{\pi_k \mathcal{N}(x_i|\mu_k, \Sigma_k)}{\sum_{j=1}^K \pi_j\mathcal{N}(x_i|\mu_j, \Sigma_j)}$$
	* 在计算$\gamma(i, k)$的时候我们假定$\mu_k$和$\Sigma_k$均已知，我们将取上一次迭代所得的值（或者初始值）
2. M步
	* $$\displaystyle
\begin{aligned}
\mu_k & = \frac{1}{N_k}\sum_{i=1}^N\gamma(i, k)x_i \\
\Sigma_k & = \frac{1}{N_k}\sum_{i=1}^N\gamma(i,
k)(x_i-\mu_k)(x_i-\mu_k)^T \\
\pi_k & = N_k/N
\end{aligned}$$

	* 其中$N_k = \sum_{i=1}^N \gamma(i, k)$，N为总数据点数。

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

### 实现

## 参考

* [漫谈 Clustering 系列](http://blog.pluskid.org/?page_id=78)<commit_msg>refresh<commit_after>--- 
layout: post 
title: Clustering杂谈
date: 2016-10-21 
categories: blog 
tags: [聚类, 机器学习] 
description: 聚类算法
--- 

# Clustering杂谈

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

### 二分k-means
 
关于k-means的变种，还有二分k-means。二分k-means的算法步骤如下：


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

## 2.k-medoids

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

$$\displaystyle\tilde{J} = \sum_{n=1}^N\sum_{k=1}^K r_{nk}\mathcal{V}(x_n,\mu_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）。

## 3.Gaussian Mixture Model（GMM）

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

$$\displaystyle
\begin{aligned}
p(x) & = \sum_{k=1}^K p(k)p(x|k) \\
& = \sum_{k=1}^K \pi_k \mathcal{N}(x|\mu_k, \Sigma_k)
\end{aligned}
$$

考虑到各个数据点的分布独立，则有：

$$p(x)=\prod_{i=1}^N p(x_i)=\displaystyle
\sum_{i=1}^N \sum_{k=1}^K \pi_k \mathcal{N}(x_i|\mu_k, \Sigma_k)
$$

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

$$\log p(x)=\displaystyle
\sum_{i=1}^N \log \left\{\sum_{k=1}^K \pi_k \mathcal{N}(x_i|\mu_k, \Sigma_k)\right\}
$$

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

1. E步
	* 对于每个数据$x_i$来说，它由第k个聚类点生成的概率为：
	$$\displaystyle
	\gamma(i, k) = \frac{\pi_k \mathcal{N}(x_i|\mu_k, \Sigma_k)}{\sum_{j=1}^K \pi_j\mathcal{N}(x_i|\mu_j, \Sigma_j)}$$
	* 在计算$\gamma(i, k)$的时候我们假定$\mu_k$和$\Sigma_k$均已知，我们将取上一次迭代所得的值（或者初始值）
2. M步
	* $$\displaystyle
\begin{aligned}
\mu_k & = \frac{1}{N_k}\sum_{i=1}^N\gamma(i, k)x_i \\
\Sigma_k & = \frac{1}{N_k}\sum_{i=1}^N\gamma(i,
k)(x_i-\mu_k)(x_i-\mu_k)^T \\
\pi_k & = \frac{N_k}{N}
\end{aligned}$$

	* 其中$N_k = \sum_{i=1}^N \gamma(i, k)$，N为总数据点数。

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

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

### 实现

```python
# -*- 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()
```

### 结果

#### 原数据

![](http://odjt9j2ec.bkt.clouddn.com/gmm-figure_1.png)

#### 聚类结果

![](http://odjt9j2ec.bkt.clouddn.com/gmm-figure_2.png)

### 数据地址

[pickle数据](https://github.com/l11x0m7/GMM/tree/master/data)

## 参考

* [漫谈 Clustering 系列](http://blog.pluskid.org/?page_id=78)
* [Regularized Gaussian Covariance Estimation](http://freemind.pluskid.org/machine-learning/regularized-gaussian-covariance-estimation/)
{'fraction_non_alphanumeric': 0.1742889012827663, 'fraction_numerical': 0.016174010039040714, 'mean_word_length': 3.445867768595041, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 4, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5973952', 'n_tokens_mistral': 6507, 'n_tokens_neox': 6857, 'n_words': 534}