本文考虑R语言的EM算法和高斯混合模型实现。

由Kaizong Ye，Coin Ge撰写

本文我们讨论期望最大化理论，应用和评估基于期望最大化的聚类。

可下载资源

数据

我们将使用mclust软件包附带的“糖尿病”数据。

data(diabetes)

summary(diabetes)

## class glucose insulin sspg
## Chemical:36 Min. : 70 Min. : 45.0 Min. : 10.0
## Normal :76 1st Qu.: 90 1st Qu.: 352.0 1st Qu.:118.0
## Overt :33 Median : 97 Median : 403.0 Median :156.0
## Mean :122 Mean : 540.8 Mean :186.1
## 3rd Qu.:112 3rd Qu.: 558.0 3rd Qu.:221.0
## Max. :353 Max. :1568.0 Max. :748.0

期望最大化（EM）

EM（Expectation-Maximization）算法是一种迭代式方法，主要应用于包含隐藏变量（latent variable）的参数估计，在无监督学习中有着广泛的应用，EM算法迭代包含两步：

利用估计的参数值来求对数似然期望（expectation）。
通过最大化对数似然期望来更新参数。

上述是EM算法的两个基本迭代部分，在实际应用中EM算法更多的看作是一种算法思想，而不是特定的算法步骤，所以接下来我将通过一些具体应用来进一步阐述EM算法的主要思想。

1、K-Means

K-Means算法是一种简单好用的聚类方法，之所以把它提到最前面来是因为我觉得同样作为无监督学习方法，K-Means在很多方面都深刻的体现着EM算法的思想。

那么首先我们来假设一个问题的形式：

input: 一组不带label的数据： $\left\{ x^{(1)},x^{(2)},...,x^{(m)} \right\}$ ，现在要对其进行聚类。

output: 聚类模型

因为与监督学习不同，数据没有带label，所以对于之前介绍的逻辑回归等分类方法就无法应用，而是需要我们从数据自身进行发掘。下面给出K-Means算法的主要步骤：

1、随机初始化聚类重心（cluster centroids）: $\mu _ { 1 } , \mu _ { 2 } , \ldots , \mu _ { k } \in \mathbb { R } ^ { n }$
2、Repeat until convergence：{
for every i, set:
$c ^ { ( i ) } : = \arg \min _ { j } \left\| x ^ { ( i ) } - \mu _ { j } \right\| ^ { 2 }$
for every j , set:
$\mu _ { j } : = \frac { \sum _ { i = 1 } ^ { m } 1 \left\{ c ^ { ( i ) } = j \right\} x ^ { ( i ) } } { \sum _ { i = 1 } ^ { m } 1 \left\{ c ^ { ( i ) } = j \right\} }$
}

可以看到，K-Means主要包含两个迭代部分，首先是最小化重心与数据距离来更新参数，接着再对属于同一类的数据求平均值来更新重心。当进行了一段时间后重心不发生变化，K-Means算法也达到了收敛。

对应的，第一步更新参数类似于EM算法中的maximization，而更新中心则对应于expectation步。虽然和真正的EM算法有着很大的差别，但是这种迭代更新的思想却是不谋而合，所以我们可以将其看作一种简化版的EM来辅助理解。

2、高斯混合模型（GMM）

接下来要介绍的高斯混合模型（Gaussian Mixture Model）是EM算法的一个重要应用，它应用广泛，而很多时候EM算法是对GMM进行参数估计的有效手段。

2.1单高斯分布

高斯分布又名正态分布，可以说是我们最常见的一种分布之一了，而且实验表明许多受多个独立因素影响的数据（如体重身高等）都满足高斯分布。

高斯分布有两个重要参数：均值 $\mu$ 以及方差 $\sigma^{2}$ ，若随机变量 $X$ 服从高斯分布，则将其写作： $X \sim N \left( \mu , \sigma ^ { 2 } \right)$ 。

一维高斯分布的概率密度函数为： $f ( x ) = \frac { 1 } { \sqrt { 2 \pi } \sigma } \exp \left( - \frac { ( x - \mu ) ^ { 2 } } { 2 \sigma ^ { 2 } } \right)$

多维高斯分布的概率密度函数为： $f _ { \mathbf { x } } \left( x _ { 1 } , \ldots , x _ { k } \right) = \frac { 1 } { \sqrt { ( 2 \pi ) ^ { k } | \mathbf { \Sigma } | } } \exp \left( - \frac { 1 } { 2 } ( \mathbf { x } - \boldsymbol { \mu } ) ^ { \mathrm { T } } \mathbf { \Sigma } ^ { - 1 } ( \mathbf { x } - \boldsymbol { \mu } ) \right)$

其中 $| \Sigma |$ 是协方差矩阵的行列式。

2.2高斯混合模型

顾名思义，高斯混合模型就是多个高斯分布线性组合在一起形成的新的分布，理论上来说，任何分布都可以由多个混合分布组合表示，但由于高斯混合模型地应用广泛，这里我们主要讨论它。

如下图所示就是两个二维高斯分布的混合：

那么如何来描述上述的分布呢？很显然使用单个高斯模型是不够准确的（从图中可以直观地看到两个不同地聚类），所以由此引入了高斯混合模型地数学表达：

$\\p ( \boldsymbol { x } ) = \sum _ { k = 1 } ^ { K } \alpha _ { k } \mathcal { N } ( \boldsymbol { x } | \boldsymbol { \mu } _ { k } , \boldsymbol { \Sigma } _ { k } ) \qquad(1)$

其中 $\alpha_{k}$ 是权重，且 $\sum _ { k = 1 } ^ { K } \alpha _ { k }=1$ ，而 $\mathcal { N } ( \boldsymbol { x } | \boldsymbol { \mu } _ { k } , \mathbf { \Sigma } _ { k } )$ 代表第k个分量模型。

如何理解GMM呢，首先我们先引入一个K维的二值变量 $z$ ，其中某个特定的元素 $z_{k}=1$ ，其余元素都为0，因此 $z_{k}$ 满足 $z _ { k } \in \{ 0,1 \}$ 且 $\sum _ { k } z _ { k } = 1$ 。

根据上述定义， $z$ 有K个不同的状态，用来表示观测数据 $x_{i}$ 属于哪个分量模型，因此每一个观测数据 $x_{i}$ 都对应一个 $z^{i}$ ，譬如 $z_{k}^{i}=1$ 表示 $x_{i}$ 来自第k个分量模型，同时变量 $z$ 的边缘概率 $p(z_{k}^{i}=1)=\alpha_{k}$ 则表示观测数据 $x_{i}$ 属于第k个分量模型的概率。而由于 $z$ 是“1 of K”的形式，因此变量 $z$ 的分布也可以写成如下形似：

$\\p ( \boldsymbol { z } ) = \prod _ { k = 1 } ^ { K }\alpha _ { k } ^ { z _ { k } }$

相应的，当对 $z$ 指定一个特定的值时，那么 $x$ 的条件概率分布就是一个特定的高斯分布：

$\\p ( \boldsymbol { x } | z _ { k } = 1 ) = \mathcal { N } ( \boldsymbol { x } | \boldsymbol { \mu } _ { k } , \mathbf { \Sigma } _ { k } )$

也可以写成：

$\\p ( \boldsymbol { x } | \boldsymbol { z } ) = \prod _ { k = 1 } ^ { K } \mathcal { N } ( \boldsymbol { x } | \boldsymbol { \mu } _ { k } , \boldsymbol { \Sigma } _ { k } ) ^ { z _ { k } }$

那么就可以得到联合概率分布 $p(x,z)=p(z)p(x|z)$ ，而变量 $x$ 的边缘概率分布就可以通过联合概率对所有可能的 $z$ 值进行求和得到，即：

$\\p ( \boldsymbol { x } ) = \sum _ { \boldsymbol { z } } p ( \boldsymbol { z } ) p ( \boldsymbol { x } | \boldsymbol { z } ) = \sum _ { k = 1 } ^ { K } \alpha _ { k } \mathcal { N } ( \boldsymbol { x } | \boldsymbol { \mu } _ { k } , \boldsymbol { \Sigma } _ { k } )\qquad(2)$

关于公式（2）第二个等式，乍一看会有些跳跃，但仔细想一想就会知道，所有可能的 $z$ 值其实就是考虑 $x$ 属于所有K个分量模型的可能性，所以体现的形式就是对所有模型加权求和。同时我们也发现经过一番推导我们得到了最初我们给出的高斯混合模型的表达式，所以我们所定义的变量 $z$ 可以看作是一个隐含变量，用于决定观测数据属于哪个模型。

2.3使用EM方法估计GMM的参数

为了方便使用EM方法估计GMM的参数，我们首先定义一个新的量 $\gamma \left( z _ { k } \right)$ ，用它来表示后验概率 $p \left( z _ { k } = 1 | \boldsymbol { x } \right)$ 。根据贝叶斯理论，我们通过先验概率 $p(\boldsymbol { z})$ 以及条件概率 $p(\boldsymbol { x }|\boldsymbol { z })$ 我们可以写出：

$\\\begin{aligned} \gamma \left( z _ { k } \right) \equiv p \left( z _ { k } = 1 | \boldsymbol { x } \right) & = \frac { p \left( z _ { k } = 1 \right) p ( \boldsymbol { x } | z _ { k } = 1 ) } { \sum _ { j = 1 } ^ { K } p \left( z _ { j } = 1 \right) p ( \boldsymbol { x } | z _ { j } = 1 ) } \\ & = \frac { \alpha _ { k } \mathcal { N } ( \boldsymbol { x } | \boldsymbol { \mu } _ { k } , \boldsymbol { \Sigma } _ { k } ) } { \sum _ { j = 1 } ^ { K } \alpha _ { j } \mathcal { N } ( \boldsymbol { x } | \boldsymbol { \mu } _ { j } , \boldsymbol { \Sigma } _ { j } ) } \end{aligned}\qquad(3)$

上面的概率公式给出的就是在已知观测样本的前提下估计隐藏变量的概率，这是我们在使用EM算法中的一个重要的过渡参数。

2.3.1极大似然函数

为了估计GMM的参数，我们要估计 $\alpha,\mu,\Sigma$ 这三个参数，所以我们首先将公式（1）的GMM模型改写成： $p ( \boldsymbol { x } | \boldsymbol { \pi } , \boldsymbol { \mu } , \boldsymbol { \Sigma } )$ 。然后由于 $x_{i}， i，2...，n$ 之间独立同分布，因此：

$\\p(X|\alpha,\mu,\Sigma)=\prod_{i=1}^{N}p(\boldsymbol { x }|\alpha,\mu,\Sigma)\qquad(4)$

为了估计参数，我们要对其进行最大化似然函数操作，为了计算方便，我们首先给出公式（4）对数似然函数：

$\\\ln p ( \boldsymbol { X } | \boldsymbol {\alpha } , \boldsymbol { \mu } , \mathbf { \Sigma } ) = \sum _ { n = 1 } ^ { N } \ln \left\{ \sum _ { k = 1 } ^ { K }\alpha _ { k } \mathcal { N } \left( \boldsymbol { x } _ { n } | \boldsymbol { \mu } _ { k } , \mathbf { \Sigma } _ { k } \right) \right\}\qquad(5)$

在最大化似然函数的前提下，我们首先来估计参数 $\mu$ ，很自然的做法就是让公式（5）对变量 $\mu_{k}$ 求导，然后令其导数为0，即：

$\\0 = \sum _ { n = 1 } ^ { K } \underbrace { \frac { \alpha _ { k } \mathcal { N } \left( \boldsymbol { x } _ { n } | \boldsymbol { \mu } _ { k } , \boldsymbol { \Sigma } _ { k } \right) } { \sum _ { j }\alpha_ { j } \mathcal { N } \left( \boldsymbol { x } _ { n } | \boldsymbol { \mu } _ { j } , \boldsymbol { \Sigma } _ { j } \right) } } _ { \gamma \left( z _ { n k } \right) } \mathbf { \Sigma } _ { k } ^ { - 1 } \left( \boldsymbol { x } _ { n } - \boldsymbol { \mu } _ { k } \right)\qquad(6)$

接着对公式（6）两边同乘 $\Sigma_{k}$ 后化简可得：

$\\\boldsymbol { \mu } _ { k } = \frac { 1 } { N _ { k } } \sum _ { n = 1 } ^ { N } \gamma \left( z _ { n k } \right) \boldsymbol { x } _ { n }\qquad(7)$

其中 $N _ { k } = \sum _ { n = 1 } ^ { N } \gamma \left( z _ { n k } \right)$ ，（化简借用公式（3））。

同样我们对 $\boldsymbol { \Sigma } _ { k }$ 求导并令其导数为0，可得：

$\\\boldsymbol { \Sigma } _ { k } = \frac { 1 } { N _ { k } } \sum _ { n = 1 } ^ { N } \gamma \left( z _ { n k } \right) \left( \mathbf { x } _ { n } - \boldsymbol { \mu } _ { k } \right) \left( \mathbf { x } _ { n } - \boldsymbol { \mu } _ { k } \right) ^ { \mathrm { T } }\qquad(8)$

最后我们要估计的参数是 $\alpha_{k}$ ，但这里需要注意的是，在一开始我们对权值 $\alpha$ 有过限制，即： $\sum _ { k = 1 } ^ { K } \alpha _ { k }=1$ ，因此在极大化似然函数时要考虑这一条件。所以我们使用拉格朗日法乘数法来求下面式子的极大值：

$\\\ln p ( \mathbf { X } | \alpha , \mu , \Sigma ) + \lambda \left( \sum _ { k = 1 } ^ { K } \alpha _ { k } - 1 \right)\qquad(9)$

然后可得：

$\\0 = \sum _ { n = 1 } ^ { N } \frac { \mathcal { N } \left( \mathbf { x } _ { n } | \boldsymbol { \mu } _ { k } , \mathbf { \Sigma } _ { k } \right) } { \sum _ { j } \alpha _ { j } \mathcal { N } \left( \mathbf { x } _ { n } | \boldsymbol { \mu } _ { j } , \mathbf { \Sigma } _ { j } \right) } + \lambda\qquad(10)$

对公式（10）前面分子分母约分后可得： $0=\sum_{n=1}^{N}{\frac{1}{\sum_{j}^{}{\alpha_{j}}}}+\lambda$ ，同时又由于 $\sum _ { k = 1 } ^ { K } \alpha _ { k }=1$ ，所以可得 $\lambda=-N$ 。接着我们对公式（10）两边同乘 $\alpha_{k}$ ，进行化简后我们可以得到：

$\\\pi _ { k } = \frac { N _ { k } } { N }\qquad(11)$

2.3.2 EM算法应用于GMM的算法步骤

根据（7）、（8）、（11）三式我们已经得到了GMM重要参数的更新迭代公式，这就是EM算法应用于高斯混合模型的实例，其过程包含E步（计算 $\gamma \left( z _ { k } \right)$ ）和M步（更新参数 $\alpha,\mu,\Sigma$ ）。接下来我们把具体的算法过程写下来：

Input：观测数据 $X$ ： $\left\{ x_{1},x_{2},...,x_{N} \right\}$
（1）初始化GMM参数：均值 $\mu_{k}$ ，协方差 $\Sigma_{k}$ ，权值 $\alpha_{k}$
（2）E步：根据当前模型的参数，计算分模型k对观测数据 $x_{j}$ 的响应度（responsibility）
$\\\gamma \left( z _ { n k } \right) = \frac { \pi _ { k } \mathcal { N } \left( \mathrm { x } _ { n } | \boldsymbol { \mu } _ { k } , \boldsymbol { \Sigma } _ { k } \right) } { \sum _ { j = 1 } ^ { K } \pi _ { j } \mathcal { N } \left( \mathbf { x } _ { n } | \boldsymbol { \mu } _ { j } , \boldsymbol { \Sigma } _ { j } \right) }$
（3）M步: 迭代更新模型的参数值：
$\\\begin{aligned} \boldsymbol { \mu } _ { k } ^ { \mathrm { new } } & = \frac { 1 } { N _ { k } } \sum _ { n = 1 } ^ { N } \gamma \left( z _ { n k } \right) \mathbf { x } _ { n } \\ \mathbf { \Sigma } _ { k } ^ { \mathrm { new } } & = \frac { 1 } { N _ { k } } \sum _ { n = 1 } ^ { N } \gamma \left( z _ { n k } \right) \left( \mathbf { x } _ { n } - \boldsymbol { \mu } _ { k } ^ { \mathrm { new } } \right) \left( \mathbf { x } _ { n } - \boldsymbol { \mu } _ { k } ^ { \mathrm { new } } \right) ^ { \mathrm { T } } \\ \pi _ { k } ^ { \mathrm { new } } & = \frac { N _ { k } } { N } \end{aligned}$
其中： $\\N _ { k } = \sum _ { n = 1 } ^ { N } \gamma \left( z _ { n k } \right)$
（4）：估算对数似然函数的值：
$\\\ln p ( \mathbf { X } | \boldsymbol { \mu } , \mathbf { \Sigma } , \boldsymbol { \pi } ) = \sum _ { n = 1 } ^ { N } \ln \left\{ \sum _ { k = 1 } ^ { K } \pi _ { k } \mathcal { N } \left( \mathbf { x } _ { n } | \boldsymbol { \mu } _ { k } , \mathbf { \Sigma } _ { k } \right) \right\}$
重复（2）、（3）步直到算法收敛（似然函数值变化或参数变化小于设定值）

期望最大化（EM）算法是用于找到最大似然的或在统计模型参数，其中该模型依赖于未观察到的潜变量最大后验（MAP）估计的迭代方法。期望最大化（EM）可能是无监督学习最常用的算法。

似然函数

似然函数找到给定数据的最佳模型。

课程

R语言数据分析挖掘必知必会

从数据获取和清理开始，有目的的进行探索性分析与可视化。让数据从生涩的资料，摇身成为有温度的故事。

立即参加 ➜

期望最大化（EM）算法

假设我们翻转硬币并得到以下内容 – 0,1,1,0,0,1,1,0,0,1。我们可以选择伯努利分布

或者，如果我们有以厘米为单位的人的身高（男性和女性）的数据。高度遵循正常的分布，但男性（平均）比女性高，因此这表明两个高斯分布的混合模型。

贝叶斯信息准则（BIC）

以糖尿病数据为例

EM集群与糖尿病数据使用mclust。

log.likelihood：这是BIC值的对数似然值

n：这是X点的数量

df：这是自由度

BIC：这是贝叶斯信息标准; 低是好的

ICL：综合完整X可能性 - BIC的分类版本。

clPairs(X,class.d)

EM的绘图命令会生成以下四个绘图：

BIC值用于选择簇的数量

聚类图

分类不确定性的图表

簇的轨道图

最受欢迎的见解

1.R语言k-Shape算法股票价格时间序列聚类

2.R语言中不同类型的聚类方法比较

3.R语言对用电负荷时间序列数据进行K-medoids聚类建模和GAM回归

4.r语言鸢尾花iris数据集的层次聚类

5.Python Monte Carlo K-Means聚类实战

6.用R进行网站评论文本挖掘聚类

7.用于NLP的Python：使用Keras的多标签文本LSTM神经网络

8.R语言对MNIST数据集分析探索手写数字分类数据

9.R语言基于Keras的小数据集深度学习图像分类

簇的散点矩阵图

可下载资源

关于作者

Kaizong Ye是拓端研究室（TRL）的研究员。在此对他对本文所作的贡献表示诚挚感谢，他在上海财经大学完成了统计学专业的硕士学位，专注人工智能领域。擅长Python.Matlab仿真、视觉处理、神经网络、数据分析。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

非常感谢您阅读本文，如需帮助请联系我们！

R语言：EM算法和高斯混合模型的实现

数据

期望最大化（EM）

1、K-Means

2、高斯混合模型（GMM）

似然函数

期望最大化（EM）算法

贝叶斯信息准则（BIC）

EM的绘图命令会生成以下四个绘图：

相关文章

R语言：EM算法和高斯混合模型的实现

数据

期望最大化（EM）

1、K-Means

2、高斯混合模型（GMM）

似然函数

期望最大化（EM）算法

贝叶斯信息准则（BIC）

EM的绘图命令会生成以下四个绘图：

相关文章

关注有关新文章的微信公众号