高斯混合聚类

简介

高斯混合聚类不同于k-means、LVQ利用原型向量刻画聚类结构，而是利用概率来刻画聚类结构。

简单来说，这种算法认为数据集中的每个样本都符合一个多元高斯分布（多元的原因是样本常是多元向量），如下

所有的样本共同符合“混合高斯分布”。混合高斯分布对应的概率密度函数是所有多元高斯分布密度函数的加权量。

多元高斯分布

若$x$服从多元高斯分布，对应概率密度函数为

$p(x) = \frac{1}{(2\pi)^{\frac{n}{2}}\lvert \Sigma \rvert^{\frac{1}{2}}}e^{-\frac{1}{2}(x-\mu)^{T}\Sigma^{-1}(x-\mu)}$，其中$x$是样本对应的向量，$\Sigma$是协方差矩阵，$\mu$是期望

为了便于理解，参照一元高斯分布

$\frac{1}{\sqrt{2\pi}\sigma}e^{\frac{-1}{2}\frac{(x - \mu)^{2}}{\sigma^{2}}}$

协方差矩阵就对应方差、多元高斯分布中期望对应一元中的期望（只是多元高斯分布中期望是一个多维向量）

所以多元高斯分布的情况，由$\Sigma$和$\mu$唯一确定

混合高斯分布

将多元高斯分布密度函数记作$p(x|\mu,\Sigma)$

可以定义混合高斯分布如下：

若$x$服从混合高斯分布，整个样本空间对应有k种多元高斯分布，对应概率密度函数$p_M = \sum^{k}{i=1} \alpha_i p(x|\mu_i,\Sigma_i)$，其中$\alpha$是混合系数，$\alpha_i$对应的实际意义是，选择第$i$个混合成分的概率，所以有$\alpha_i > 0$且$\sum^{k}{i=1} \alpha_i = 1$；而$p(x|\mu_i,\Sigma_i)$对应第i个混合成分的概率密度

样本属于某混合成分的概率

令数据集$D = {x_1,x_2,…,x_m}$随机变量$z_j \in {1,2,…,k}$，$z_j$表征样本$x_j$属于哪个混合成分。

对于$x_j$并不确定的情况下，$z_j$的先验分布为

$p(z_j = i) = \alpha_i, i = {1,2,…,k}$

根据贝叶斯定理，当$x_j$确定时，$z_j$的后验分布记作$p_M(z_j = i|x_j)$，为

$p_M(z_j = i|x_j) = \frac{p(z_j = i)p_M(x_j|z_j = i)}{p_M(x_j)}$
其中，$p(z_j = i) = \alpha_i$；$p_M(x_j|z_j = i) = p(x|\mu_i,\Sigma_i)$，因为当确定$z_j = i$时除了$\alpha_i = 1$其它$\alpha_1,…\alpha_i-1,\alpha_i+1,…,\alpha_k$均为$0$；而$p_M(x_j)$即混合高斯分布的概率密度函数

将$p_M(z_j = i|x_j)$简记作$\gamma_{ji}$，这个概率就是$x_j$的分布为$\mu_i,\Sigma_i$所对应的多元高斯分布的概率。

于是我们可以找到令$\gamma_{ji}$最大的$i \in {1,2,…,k}$，令$\lambda_j = argmax_{i \in {1,2,…,k}} \gamma_{ji}$，$\lambda_j$即$x_j$的标签。对于每个都进行相同的操作，整个数据集便可被划分为k个多元高斯分布对应的k个簇。

模型训练

目的

对于k-means、LVQ的模型训练实际上就是通过训练集获得对应的原型向量，有了原型向量便有了划分为簇的依据，也就完成了模型的训练。

而对于高斯混合聚类，根据前面的描述，我们划分为簇的重要依据就是$\gamma_{ji}$，进一步说实际上是计算$\gamma_{ji}$的依据，根据计算公式可知，这个依据实际上就是决定高斯混合分布的参数$\mu_i,\Sigma_i,\alpha_i, i \in {1,2,…,k}$。

于是训练的目的实际上就是要得到它们的值。

过程

对于模型参数${(\alpha_i,\mu_i,\Sigma_i)|1 \le i \le k}$，我们采用极大似然的方法进行求解。

这里引用南瓜书中的一句话：
“对于每个样本$x_j$来说，它出现的概率是$p_M(x_j)$既然现在训练集D中确实出现了$x_j$，我们当然希望待求解的参数${(\alpha_i,\mu_i,\Sigma_i)|1 \le i \le k}$能够使这种可能性$p_M(x_j)$最大”

于是根据极大似然方法，我们令
$LL(D) = ln(\prod^{m}{j=1} p_M(x_j)) = \sum^{m}{j=1} ln(\sum^k_{i=1} \alpha_i \cdot p(x_j|\mu_i,\Sigma_i))$
为对数似然函数，将其最大化得到对应的参数，就是我们要求解的模型参数。

结果

经过一系列数学运算，我们可以得得到如下结果：

$\mu_i = \frac{\sum^{m}_{j=1} \gamma_{ji} x_j}{\sum^{m}_{j=1} \gamma_{ji}}$

$\Sigma_i = \frac{\sum^{m}_{j=1} \gamma_{ji} (x_j - \mu_i)(x_j - \mu_i)^{T}}{\sum^{m}_{j=1}} \gamma_{ji}$

$\alpha_i = \frac{1}{m} \sum^{m}_{j=1} \gamma{ji}$

$i = 1,2,…,k$

注意结果当中，每一个参数的计算都要用到$\gamma{ji}$，而问题在于计算$\gamma{ji}$的时候又要用到参数本身，这似乎是一个循环的无从下手的问题。这种情况下，我们利用EM算法来求解。

EM算法

EM算法(Expectation-Maximization)称为“期望最大化算法”，这种算法最开始应用于：使用极大似然法对模型参数进行估计，但是已知的样本中存在还没有“观测”的变量，这种变量称为隐变量，它的值是不确定的。

令$X,Z,\Theta$分别为已观测变量集、隐变量集、参数集，则应最大化对数似然$LL(\Theta|X,Z) = lnP(X,Z|\Theta)$，但是Z是隐变量，所以无法直接求解。

EM算法可以用于估计隐变量，并在这个过程中对参数做最大似然估计。

其基本的思想是这样的：

1. 初始化参数$\Theta$，根据参数去估计隐变量的概率分布，并利用此概率分布求得隐变量的期望——E步
2. 将隐变量的期望作为我们观测到的隐变量本身，于是此时所有样本都已被观测，对$\Theta$做极大似然估计——M步

不断重复上述两个过程——迭代，直到满足退出条件，例如$\Theta$收敛

贴一篇介绍EM算法的博客：
CSDN

算法流程

结合EM算法，整个高斯混合聚类算法流程如下：

输入：

1. 样本集$D = {x_1,x_2,…,x_m}$
2. 高斯混合成分的个数k（当然也就是希望划分出的簇的个数）

过程：

1. 初始化高斯混合分布的参数${(\alpha_i,\mu_i,\Sigma_i)|i = 1,2…,k}$
2. repeat：
3. 对每个样本$x_j,j = 1,2,…,m$估计其属于第$i,i = 1,2,…k$个成分的概率：$\gamma_{ji}$
4. 利用公式
   $\mu_i = \frac{\sum^{m}_{j=1} \gamma_{ji} x_j}{\sum^{m}_{j=1} \gamma_{ji}}$，
   $\Sigma_i = \frac{\sum^{m}_{j=1} \gamma_{ji} (x_j - \mu_i)(x_j - \mu_i)^{T}}{\sum^{m}_{j=1}} \gamma_{ji}$，
   $\alpha_i = \frac{1}{m} \sum^{m}_{j=1} \gamma{ji}$，$i = 1,2,…,k$对参数进行更新
5. until:收敛条件（达到一定轮数 or 参数收敛）
6. 求解$x_j,j = 1,2,…,m$的标记$\lambda_j = agrmax_{i = 1,2,…,k} \gamma_{ji}$
7. 根据标记划分为对应的簇

代码实现

Data.py:数据集部分

1
2
3
4
5
6
7
8
9

D = [
    [0, 0],
    [0.697, 0.460], [0.774, 0.376], [0.634, 0.264], [0.608, 0.318], [0.556, 0.215],
    [0.403, 0.237], [0.481, 0.149], [0.437, 0.211], [0.666, 0.091], [0.243, 0.267],
    [0.245, 0.057], [0.343, 0.099], [0.639, 0.161], [0.657, 0.198], [0.360, 0.370],
    [0.593, 0.042], [0.719, 0.103], [0.359, 0.188], [0.339, 0.241], [0.282, 0.257],
    [0.748, 0.232], [0.714, 0.346], [0.483, 0.312], [0.478, 0.437], [0.525, 0.369],
    [0.751, 0.489], [0.532, 0.472], [0.473, 0.376], [0.725, 0.445], [0.446, 0.459]
]  # 数据集，1~30，0索引不使用

main.py:主函数部分

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

import Gauss
import Data


def main():
    print("高斯混合聚类 的结果", end="\n")
    res = Gauss.gauss(Data.D)
    for i in range(3):
        print(len(res[i]), end="\n")
        print(res[i])


if __name__ == "__main__":
    main()

Guass.py:高斯混合聚类算法部分

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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106

import numpy as np


def multi_gauss_distri_p(sigma, mu, n, x):
    # 计算多元高斯分布下取得x的概率，n是维度
    vec_mu = np.array(mu)
    vec_x = np.array(x)
    t1 = vec_x - vec_mu
    #    t1 = t1.T
    det_sigma = np.linalg.det(sigma)
    if(det_sigma < 1e-10):  # 当行列式过小时，添加一个较小的正则化项
        sigma += np.eye(sigma.shape[0]) * 1e-6  # 添加正则化，避免奇异矩阵
    t2 = np.linalg.inv(sigma)
    t3 = vec_x - vec_mu
    t3 = t3.T  # 西瓜书上的公式里没有转置的向量默认是列向量
    log_p = -0.5 * (np.dot(np.dot(t1, t2), t3) + np.linalg.slogdet(sigma)[1] + n * np.log(2 * np.pi))
    p = np.exp(log_p)
#    p = 1 / ((2 * np.pi) ** (n / 2) * np.linalg.det(sigma) ** (1 / 2)) * np.e ** (-0.5 * np.dot(np.dot(t1, t2), t3)
#    print("p:", p, end="\n")
    if p < 1e-10:
        p = 1e-6  # 避免数值问题
    return p


def new_mu_i(D, lamda_, i):
    up_sum = [0, 0]
    down_sum = 0
    for j in range(1, 31):
        up_sum[0] += lamda_[j][i] * D[j][0]
        up_sum[1] += lamda_[j][i] * D[j][1]
        down_sum += lamda_[j][i]
    new_mu = [up_sum[0] / down_sum, up_sum[1] / down_sum]
    return new_mu


def new_sigma_i(D, mu, lamda_, i):
    vec_x = np.array(D[i])
    vec_mu = np.array(mu)
    up_sum = np.array([[0.0, 0.0], [0.0, 0.0]])
    down_sum = 0
    for j in range(1, 31):
        t = vec_x - vec_mu
        up_sum += lamda_[j][i] * (np.outer(t.T, t))
        down_sum += lamda_[j][i]
    new_sigma = up_sum / down_sum
    return new_sigma


def new_alpha_i(lambda_, i):
    up_sum = 0
    for j in range(1, 31):
        up_sum += lambda_[j][i]
    new_alpha = up_sum / 30
    return new_alpha


def gauss(D: [[float]]):
    #  定义混合高斯分布参数
    sigma = []  # 协方差矩阵
    alpha = []  # 混合权重
    mu = []  # 期望
    #  初始化参数
    for i in range(0, 3):
        sigma.append(np.array([[0.1, 0.0], [0.0, 0.1]]))
        alpha.append(1/3)
    for i in range(0, 3):
        if i == 0:
            mu.append(D[6])
        elif i == 1:
            mu.append(D[22])
        elif i == 2:
            mu.append(D[27])
    #  定义xj属于第i种多元高斯分布的概率
    lambda_ = []
    for j in range(0, 31):
        lambda_.append([[], [], []])
    for _ in range(0, 100):  # 迭代
        for j in range(1, 31):  # 依次对每个xj计算lambda_ji,i = 0，1，2，对应三种多元高斯分布
            for i in range(0, 3):
                t_sum = 0
                for l in range(0, 3):
                    t_sum += alpha[l] * multi_gauss_distri_p(sigma[l], mu[l], 2, D[j])
                lambda_ji = alpha[i] * multi_gauss_distri_p(sigma[i], mu[i], 2, D[j]) / t_sum
                lambda_[j][i] = lambda_ji
        for i in range(0, 3):  # 更新混合高斯分布参数
            mu[i] = new_mu_i(D, lambda_, i)
            sigma[i] = new_sigma_i(D, mu[i], lambda_, i)
            alpha[i] = new_alpha_i(lambda_, i)

    result = [[], [], []]  # 最终的划分结果
    for j in range(1, 31):  # 依次对每个xj计算lambda_ji,i = 0，1，2，对应三种多元高斯分布
        for i in range(0, 3):
            t_sum = 0
            for l in range(0, 3):
                t_sum += alpha[l] * multi_gauss_distri_p(sigma[l], mu[l], 2, D[j])
            lambda_ji = alpha[i] * multi_gauss_distri_p(sigma[i], mu[i], 2, D[j]) / t_sum
            lambda_[j][i] = lambda_ji

        flag = 0
        ll = lambda_[j][0]
        for i in range(1, 3):
            if ll < lambda_[j][i]:
                flag = i
        result[flag].append(D[j])
    return result

进行这一部分的时候让我最大的感悟是，在使用计算机进行数据处理的时候，很可能会出现数据的溢出问题，非常大、非常小的数都在计算机的表示范围之外，就会带来问题。例如在进行协方差矩阵更新的时候，有时它的行列式值虽然不是0，但是已经非常小了，计算机会默认其为0。同样，可能作为除数的数也是一样的，如果太小变为0就会发送除以0的错误。
这称为数值稳定性问题

因此代码中有一些涉及处理这些问题的地方（主要是在求第i个多元高斯分布中取得xj这个值的概率的时候，一处是正则化、一处是对数化并添加过小的判断），尽管我现在并不清楚这样处理是否有问题。但是抛开这些细节，作为一次上手的练习，整个混合高斯聚类算法是正确的。
不过值得一提的是，EM算法的收敛与迭代的次数没有必然的关系，通常应该使用参数变化的情况来决定是否结束迭代

最后的结果如下：

高斯混合聚类 的结果
2
[[0.608, 0.318], [0.359, 0.188]]
0
[]
28
[[0.697, 0.46], [0.774, 0.376], [0.634, 0.264], [0.556, 0.215], [0.403, 0.237], [0.481, 0.149], [0.437, 0.211], [0.666, 0.091], [0.243, 0.267], [0.245, 0.057], [0.343, 0.099], [0.639, 0.161], [0.657, 0.198], [0.36, 0.37], [0.593, 0.042], [0.719, 0.103], [0.339, 0.241], [0.282, 0.257], [0.748, 0.232], [0.714, 0.346], [0.483, 0.312], [0.478, 0.437], [0.525, 0.369], [0.751, 0.489], [0.532, 0.472], [0.473, 0.376], [0.725, 0.445], [0.446, 0.459]]

经过我的观察，第3次迭代之后划分结果基本上就稳定了，第1次的时候分得比较均匀。这或许是数据的原因。（也有可能是我对于数值稳定性的处理不好）

如果将正则化项改大一些结果如下：

高斯混合聚类 的结果
12
[[0.608, 0.318], [0.556, 0.215], [0.403, 0.237], [0.481, 0.149], [0.437, 0.211], [0.243, 0.267], [0.245, 0.057], [0.343, 0.099], [0.359, 0.188], [0.339, 0.241], [0.282, 0.257], [0.483, 0.312]]
10
[[0.634, 0.264], [0.666, 0.091], [0.639, 0.161], [0.657, 0.198], [0.593, 0.042], [0.719, 0.103], [0.748, 0.232], [0.714, 0.346], [0.751, 0.489], [0.725, 0.445]]
8
[[0.697, 0.46], [0.774, 0.376], [0.36, 0.37], [0.478, 0.437], [0.525, 0.369], [0.532, 0.472], [0.473, 0.376], [0.446, 0.459]]