EM算法求解GMM的python实现——基于kMeans实现参数初始化

EM算法求解GMM的python实现——基于kMeans实现参数初始化

本文重点参考该篇博文：
版权声明：本文为CSDN博主「deephub」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/deephub/article/details/129791864

一、理论知识

1.1、高斯混合模型（GMM）及期望最大（EM）算法

1.1.1、GMM

（1）基本概念

• 高斯混合模型（Gaussian Mixture Model, GMM）是单一高斯概率密度函数的延伸，GMM能够平滑地近似任意形状的密度分布。

• 高斯混合模型种类有单高斯模型（Single Gaussian Model, SGM）和高斯混合模型（Gaussian Mixture Model, GMM）两类。

• 根据高斯概率密度函数（Probability Density Function, PDF）参数不同，每一个高斯模型可以看作一种类别，输入一个样本x，即可通过PDF计算其值，然后通过一个阈值来判断该样本是否属于高斯模型。

• SGM适合于仅有两类别问题的划分，而GMM由于具有多个模型，划分更为精细，适用于多类别的划分，可以应用于复杂对象建模。

  • k-means聚类模型非常简单并且易于理解，但是实际应用中，k-means的非概率性和它仅根据到簇中心点的距离来指派将导致性能低下；
  • 高斯混合模型可以看作是KMeans的一个扩展（“概率版本的KMeans”），是一种非常强大的聚类评估工具。

（2）模型参数估计

详细内容可阅读该文章：

• 高斯混合模型（GMM） - 戴文亮的文章 - 知乎 https://zhuanlan.zhihu.com/p/30483076
• https://blog.csdn.net/jojozhangju/article/details/19182013

首先定义如下信息：

• $x_j$表示第$j$ 个观测数据，$j=1,2,...,N$ ；
• $k$表示混合模型中第$k$个子模型（第$k$类），$k=1,2,...,K$ ；
• ${\alpha }_k$表示观测数据属于第$k$个子模型的概率，${\alpha }_k\ge 0,{\sum }_{k=1}^K{\alpha }_k=1$ ；
• $\varphi (x|{\theta }_k)$是第$k$个子模型的高斯分布密度函数，${\theta }_k=({\mu }_k,{\delta }_k^2)$ ，其展开形式就是正态分布的概率密度函数；
• ${\gamma }_{jk}$表示第$j$个观测数据属于第$k$个子模型的概率。

高斯混合模型的概率分布为：
$$P(x|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (x|{\theta }_k)$$
对于这个模型而言，参数${\theta }_k=({\mu }_k,{\delta }_k,{\alpha }_k)$，也就是每个子模型的期望、方差（或协方差）、在混合模型中发生的概率。

（2.1）对于SGM（样本类别已知）或者样本类别已知的GMM

• 可以直接求解模型参数${\mu }_k,{\delta }_k,{\alpha }_k$

设属于第$k$个子模型的样本集合为$S(k)$，样本数量为$N_k$
$${\alpha }_k=\frac{N_k}{N}$$

$${\mu }_k=\frac{1}{N_k}\sum _{x\in S(k)}x$$

$$\sum _k=\frac{1}{N_k}\sum _{x\in S(k)}(x-{\mu }_k)(x-{\mu }_k{)}^T,或者表示为{\delta }_k=\frac{1}{N_k}\sum _{x\in S(k)}(x-{\mu }_k)(x-{\mu }_k{)}^2$$

（2.2）对于样本类别未知的GMM

目标：希望找到一组参数$\theta$，使得生成已有的$N$个样本的概率最大，这个概率就是$L(\theta )={\prod }_{j=1}^{N}P(x_j|\theta )$

由于每个点发生的概率都很小，乘积会变得极其小，不利于计算和观察，因此通常用最大似然估计（Maximum Log-Likelihood）来计算

• Log函数具备单调性，不会改变极值的位置，同时在 0-1 之间输入值很小的变化可以引起输出值相对较大的变动

• 所求概率由$L(\theta )={\prod }_{j=1}^{N}P(x_j|\theta )$变为$logL(\theta )={\sum }_{j=1}^{N}logP(x_j|\theta )$

• 对于高斯混合模型，$logL(\theta )={\sum }_{j=1}^{N}logP(x_j|\theta )={\sum }_{j=1}^{N}log({\sum }_{k=1}^K{\alpha }_{k}P(x_j|{\theta }_k))$

观察该公式可知，无法直接求导确定使 likelihood 最大的参数

• 对于每个观测数据点，事先并不知道它是属于哪个子分布的（hidden variable），对应$log$里面的求和项${\sum }_{k=1}^K{\alpha }_{k}P(x_j|{\theta }_k)$
• 对于每个子模型都有未知的${\mu }_k,{\delta }_k,{\alpha }_k$，直接求导无法计算
• 需要通过迭代的方法求解（EM算法）。

1.1.2、EM算法

EM 算法是一种迭代算法，1977 年由 Dempster 等人总结提出，用于含有隐变量（Hidden variable）的概率模型参数的最大似然估计。

每次迭代包含两个步骤：

1. E-step：求期望$E({\gamma }_{jk}|X,\theta )$for all $j=1,2,...,N$
2. M-step：求极大，计算新一轮迭代的模型参数${\mu }_k,{\delta }_k,{\alpha }_k$for all $k=1,2,...,K$

一般性的 EM 算法：通过 Jensen 不等式得出似然函数的下界 Lower bound，通过极大化下界做到极大化似然函数

相关内容可参考链接：EM算法存在的意义是什么？ - Mark的回答 - 知乎 https://www.zhihu.com/question/40797593/answer/3058344957

在高斯混合模型里应用EM算法估计模型参数

1. 初始化参数${\mu }_k,{\delta }_k,{\alpha }_k$for all $k=1,2,...,K$

2. E-step：依据当前参数，计算每个数据$x_j$来自子模型$k$的可能性${\gamma }_{jk}$

   • 计算各个样本属于各个类别的概率（根据“数据分布”进行“样本分类”）

   • 公式推导：可参见1.2贝叶斯公式章节

     • $$P(B_K|A)=\frac{P(A|B_k)P(B_k)}{P(A)}=\frac{P(A|B_k)P(B_k)}{{\sum }_{i=1}^{n}P(A|B_i)P(B_i)}$$

   • $${\gamma }_{jk}=\frac{{\alpha }_k\varphi (x_j|{\theta }_k)}{{\sum }_{k=1}^K{\alpha }_k\varphi (x_j|{\theta }_k)},j=1,2,...,N;k=1,2,...,K$$

     • 其中，${\sum }_{k=1}^K{\alpha }_k\varphi (x_j|{\theta }_k)$表示根据各个子模型的当前参数出现数据$x_j$的概率和

3. M-step：计算模型参数，用于新一轮迭代

   • 计算模型参数（根据“样本分类”拟合“数据分布”）

   • for all $k=1,2,...,K$

   • $${\mu }_k=\frac{{\sum }_{j=1}^N{\gamma }_{jk}{x}_j}{{\sum }_{j=1}^N{\gamma }_{jk}}$$

   • $$\sum _k=\frac{{\sum }_{j=1}^N{\gamma }_{jk}(x-{\mu }_k)(x-{\mu }_k{)}^T}{{\sum }_{j=1}^N{\gamma }_{jk}},或者表示为{\delta }_k=\frac{{\sum }_{j=1}^N{\gamma }_{jk}(x-{\mu }_k)(x-{\mu }_k{)}^2}{{\sum }_{j=1}^N{\gamma }_{jk}}$$

   • $${\alpha }_k=\frac{{\sum }_{j=1}^N{\gamma }_{jk}}{N}$$

4. 重复计算 E-step 和 M-step 直至收敛：$||{\theta }_{i+1}-{\theta }_i||<ϵ$，$ϵ$是一个很小的正数

   • 经过一次迭代之后参数变化非常小，则停止迭代（也可以直接设置迭代次数作为收敛条件）

其中算法伪代码如下图：
图片来源：https://blog.csdn.net/weixin_45488228/article/details/102463264

EM算法估计GMM模型参数的示例：

【[5分钟学算法] #06 EM算法 你到底是哪个班级的】 https://www.bilibili.com/video/BV1RT411G7jJ/?share_source=copy_web&vd_source=c71860a35ae6f457f141a44ee4b38133

1.2、贝叶斯公式

1.2.1、乘法公式

对于任意两个事件$A$、$B$，有
$$P(AB)=P(A|B)P(B)=P(B|A)P(A)$$
证明：

由条件概率的定义可证，即任意两个事件同时发生的概率等于一个事件发生的概率乘以该事件已发生条件下另一事件发生的概率，要求$P(A)>0,P(B)>0$

1.2.2、全概率公式

设$B_1,B_2,...,B_n$是基本空间$\Omega$的一个分割，则对$\Omega$中任一事件$A$，有
$$P(A)=\sum _{i=1}^{n}P(A|B_i)P(B_i)$$
证明：

已知，$A=A\Omega =A(\bigcup _{i=1}^n{B}_i)=\bigcup _{i=1}^{n}A{B}_i$，由$B_1,B_2,...,B_n$互不相容，可推导出$A{B}_1,A{B}_2,...,A{B}_n$也互不相容，因此
$$P(A)=P(\bigcup _{i=1}^{n}A{B}_i)=\sum _{i=1}^{n}P(A{B}_i)=\sum _{i=1}^{n}P(A|B_i)P(B_i)$$

1.2.3、贝叶斯公式

设事件$B_1,B_2,...,B_n$是基本空间$\Omega$的一个分割，发生的概率皆已知且为正，设$A$为$\Omega$中一个事件，发生概率为正，且在诸$B_i$给定下事件$A$的条件概率$P(A|B_1),P(A|B_2)...P(A|B_n)$可通过实验等手段获得，则在$A$给定下，事件$B_k$的条件概率为
$$P(B_K|A)=\frac{P(A|B_k)P(B_k)}{P(A)}=\frac{P(A|B_k)P(B_k)}{{\sum }_{i=1}^{n}P(A|B_i)P(B_i)}$$

二、代码实现

2.1、E-step

def expectation_step(class_num, samples, possibility, means, variances):
    """
    E-step，计算各个样本属于各个类别的概率（根据“数据分布”进行“样本分类”）
    :param class_num: 类别数
    :param samples: 样本
    :param possibility: 权重（样本属于各个类别的概率）
    :param means: 均值
    :param variances: 方差
    :return: 各个样本属于各个类别的概率
    """
    # 概率初始化
    weights = np.zeros((class_num, len(samples)))
    # 根据各个类别的特征值求概率密度
    for class_index in range(class_num):
        # 求概率密度函数
        weights[class_index] = norm(loc=means[class_index], scale=np.sqrt(variances[class_index])).pdf(samples)

    # 计算不同类别出现某个样本的概率和（贝叶斯公式里的“分母项”）
    weight_sum = np.sum([weights[class_index] * possibility[class_index] for class_index in range(class_num)], axis=0)

    for class_index in range(class_num):
        # 计算每个样本属于每种类别的概率，基于贝叶斯公式
        weights[class_index] = weights[class_index] * possibility[class_index] / weight_sum

    return weights

2.2、M-step

def maximization_step(class_num, samples, weights, possibility, means, variances):
    """
    M-step，计算权重以及各个类别的特征值（根据“样本分类”拟合“数据分布”）
    :param class_num: 类别数
    :param samples: 样本
    :param weights: 各个样本属于各个类别的概率
    :param possibility: 权重（样本属于各个类别的概率）
    :param means: 均值
    :param variances: 方差
    :return: 权重、均值、方差
    """
    for class_index in range(class_num):
        # 计算每类的均值mu
        means[class_index] = np.sum(weights[class_index] * samples) / np.sum(weights[class_index])

        # 计算每类的方差sigma
        variances[class_index] = np.sum(weights[class_index] * (samples - means[class_index]) ** 2) / np.sum(weights[class_index])

        # 计算样本属于各个类别的权重
        possibility[class_index] = np.mean(weights[class_index])

    return possibility, means, variances

2.3、使用KMeans进行参数初始化

def random_init(samples, class_num, init_method='kmeans'):
    """
    初始化权重（样本属于各个类别的概率）、各个分布的均值、方差
    init_method='kmeans'，使用kmeans进行样本初始分类，确定权重、均值、方差
    init_method='random'，使用均等权重，随机确定均值、方差
    :param samples: 样本
    :param class_num: 类别数
    :return: 权重、均值、方差
    """
    if init_method == 'random':
        # 初始化样本属于各个类别的权重（均相等）
        possibility = np.ones(class_num) / class_num
        # 初始化各个分布的均值、方差
        means = np.random.choice(samples, class_num)
        variances = np.random.random_sample(class_num)

        return possibility, means, variances

    # 各样本到各质心的距离
    distance_matrix = np.zeros((class_num, len(samples)))
    # 随机选择质心
    center = np.random.choice(samples, class_num)
    # 迭代
    total_epoch = 100
    for epoch in range(total_epoch):
        # 计算距离（一维数据直接以差值的绝对值作为距离）
        for class_index in range(class_num):
            distance_matrix[class_index] = np.abs(samples - center[class_index])

        # 确定各个样本所属类别
        sample_class = np.argmin(distance_matrix, axis=0)
        # 更新质心
        for class_index in range(class_num):
            center[class_index] = np.mean([samples[i] for i, c in enumerate(sample_class) if c == class_index])

    # 根据各个类别的样本数据，计算possibility, means, variances
    possibility, means, variances = [], [], []
    # 计算距离（一维数据直接以差值的绝对值作为距离）
    for class_index in range(class_num):
        distance_matrix[class_index] = np.abs(samples - center[class_index])

    # 确定各个样本所属类别
    sample_class = np.argmin(distance_matrix, axis=0)
    for class_index in range(class_num):
        possibility.append(len([i for i in sample_class if i == class_index]) / len(samples))
        means.append(np.mean([samples[i] for i, c in enumerate(sample_class) if c == class_index]))
        variances.append(np.mean([(samples[i] - means[class_index]) ** 2 for i, c in enumerate(sample_class) if c == class_index]))

    return possibility, means, variances

sklearn提供的KMeans非常完备，实际使用时应优先使用（相比于个人手写的算法）。

2.4、使用scikit-learn提供的GMM

from sklearn.mixture import GaussianMixture

# 使用已有的库函数
gmm = GaussianMixture(n_components=3)
gmm.fit(samples.reshape(-1, 1))
# 输出各个分布的均值
print("Means:")
print(gmm.means_)

# 输出各个分布的协方差矩阵
print("Covariances:")
print(gmm.covariances_)

# 输出各个分布的权重
print("Weights:")
print(gmm.weights_)

初始化GaussianMixture()时，常用参数如下：

• n_components：类别数（默认为1）；
• tol：模型的停止标准，当下限平均增益低于tol参数时（默认为0.001），EM迭代将停止；
• init_params：用于初始化权重的方法，随机初始化（‘random’）或者使用聚类算法初始化（‘kmeans’，默认方法）；

sklearn提供的GMM非常完备，实际使用时应优先使用（相比于个人手写的算法）。

三、参考链接

代码文件：

1. https://download.csdn.net/download/weixin_42639395/88136996?spm=1001.2014.3001.5503

GMM相关：

1. https://blog.csdn.net/deephub/article/details/129791864
2. 高斯混合模型（GMM） - 戴文亮的文章 - 知乎 https://zhuanlan.zhihu.com/p/30483076
3. https://blog.csdn.net/jojozhangju/article/details/19182013

EM算法相关：

1. EM算法存在的意义是什么？ - Mark的回答 - 知乎 https://www.zhihu.com/question/40797593/answer/3058344957
2. https://blog.csdn.net/weixin_45488228/article/details/102463264
3. 【[5分钟学算法] #06 EM算法 你到底是哪个班级的】 https://www.bilibili.com/video/BV1RT411G7jJ/?share_source=copy_web&vd_source=c71860a35ae6f457f141a44ee4b38133