【机器学习】EM算法在高斯混合模型学习中的应用

前言

EM算法，此博客介绍了 $E M$ 算法相关理论知识，看本篇博客前先熟悉 $E M$ 算法。

本篇博客打算先从单个高斯分布说起，然后推广到多个高斯混合起来，最后给出高斯混合模型参数求解过程。

单个高斯分布

假如我们有一些数据，这些数据来自同一个高斯分布（独立同分布），那个我们如何根据这些数据估计出这个高斯分布的参数呢？我们知道只要知道高斯分布的参数 $θ = {μ, σ^{2}}$ 就能确定此高斯分布。

从上图中，我们要想知道数据是来自哪个高斯分布，我们就要知道高斯分布的参数，直观上可以认定数据来自参数为

θ_{1}

的高斯分布，然而这毕竟是我们直观上的，那么我们应该如何根据数据估计高斯分布的参数呢？

假设数据为 $X = {x_{1}, x_{2}, . . ., x_{N}} ， x_{i}$ 独立同分布 $p (X | θ)$ ，其中 $θ = {μ, σ^{2}}$ ;

由贝叶斯公式知：

p (θ | X) \propto p (X | θ) p (θ)

p (θ | X)

为后验概率，

p (X | θ)

为似然度，

p (θ)

为先验概率。

关于贝叶斯公式的，可参考我之前的博客，里面有提到：链接

不加上先验概率，就是极大似然估计；
加上先验概率，就是极大后验概率估计
这里我们只介绍极大似然估计，极大后验概率估计类似

（1）写出对数似然函数

L (θ | X) = l o g [p (X | θ)] = \sum_{i = 1}^{N} l o g p (x_{i} | θ) = \sum_{i = 1}^{N} l o g [\frac{1}{\sqrt{2 π} σ} e x p (- \frac{(x_{i} - μ)^{2}}{2 σ^{2}})]

（2）求极大似然估计

分别对 $μ, σ$ 求偏导，并令为0：

\frac{\partial L (θ | X)}{\partial μ} = \frac{\partial (\sum_{i = 1}^{N} l o g [\frac{1}{\sqrt{2 π} σ} e x p (- \frac{(x_{i} - μ)^{2}}{2 σ^{2}})])}{\partial μ} = \frac{\partial (\sum_{i = 1}^{N} l o g [e x p (- \frac{(x_{i} - μ)^{2}}{2 σ^{2}})])}{\partial μ} = \frac{\partial (\sum_{i = 1}^{N} - \frac{(x_{i} - μ)^{2}}{2 σ^{2}})}{\partial μ} = - \sum_{i = 1}^{N} \frac{(x_{i} - μ)}{σ^{2}} = 0

由此得到

μ

的似然估计为：

μ_{M L E} = \frac{1}{N} \sum_{i = 1}^{N} x_{i}

\frac{\partial L (θ | X)}{\partial σ} = \frac{\partial (\sum_{i = 1}^{N} l o g [\frac{1}{\sqrt{2 π} σ} e x p (- \frac{(x_{i} - μ_{M L E})^{2}}{2 σ^{2}})])}{\partial σ} = \frac{\partial \sum_{i = 1}^{N} l o g \frac{1}{\sqrt{2 π} σ}}{\partial σ} + \frac{\partial (\sum_{i = 1}^{N} l o g [e x p (- \frac{(x_{i} - μ_{M L E})^{2}}{2 σ^{2}})])}{\partial σ} = \sum_{i = 1}^{N} - (\frac{1}{\sqrt{2 π} σ}) \sqrt{2 π} + \frac{\partial (\sum_{i = 1}^{N} - \frac{(x_{i} - μ_{M L E})^{2}}{2 σ^{2}})}{\partial σ} = - N + \frac{\sum_{i = 1}^{N} (x_{i} - μ_{M L E})^{2}}{σ^{2}} = 0

其中用到

\frac{d l o g (\frac{1}{x})}{d x} = - \frac{1}{x}

得到

σ^{2}

的似然估计为：

σ^{2} = \frac{\sum_{i = 1}^{N} (x_{i} - μ_{M L E})^{2}}{N}

综上我们可以表述为：

θ = a r g max_{θ} [\sum_{i = 1}^{N} l o g N (x_{i} | μ, σ^{2})]

μ \to \frac{\partial L (μ, σ^{2} | X)}{\partial μ} σ^{2} \to \frac{\partial L (μ, σ^{2} | X)}{\partial σ}

以上问题是只有一个高斯分布的，如果不止一个高斯分布呢？

高斯混合模型

同样，我们看下图：

上图是一个高斯混合模型，此处只画了两个高斯分布，可以是多个高斯分布。

如果我们知道每一个数据属于哪一个高斯分布，就会很容易求解，但是我们不可能都知道的，这时随便一个数据点，我们应该如何判断它是哪个高斯分布产生的呢？

这里我们引入 $α_{k}$ 表示属于第 $k$ 个高斯分布的权重，并满足 $\sum_{k = 1}^{K} α_{k} = 1$

这样我们可以得到：

p (X | θ) = \sum_{k = 1}^{K} α_{k} N (μ_{k}, σ_{k}^{2}) \sum_{i = k}^{K} α_{k} = 1

其中 $N (μ_{k}, σ_{k}^{2})$ 是高斯分布， $α_{k}$ 是系数

给出定义：

高斯混合模型是指具有如下形式的概率分布模型：

p (x | θ) = \sum_{k = 1}^{K} α_{k} ϕ (x | θ_{k})

其中，

α_{k}

是系数，

α_{k} ⩾ 0 \sum_{k = 1}^{K} α_{k} = 1

；

ϕ (x | θ_{k})

是高斯分布密度，

θ_{k} = (μ_{k}, σ_{k}^{2})

，

ϕ (x | θ_{k}) = \frac{1}{\sqrt{2 π} σ_{k}} e x p (- \frac{(x - μ_{k})^{2}}{2 σ_{k}^{2}})

称为第

k

个分模型。

高斯混合模型参数估计的 $E M$ 算法

假设观测数据 $x_{1}, x_{2}, . . ., x_{N}$ 由高斯混合模型生成，

p (x | θ) = \sum_{k = 1}^{K} α_{k} ϕ (x | θ_{k})

其中，

θ = (α_{1}, α_{2}, . ., α_{k}; θ_{1}, θ_{2}, . ., θ_{k})

，我们用

E M

算法估计高斯混合模型的参数

θ

1. 明确隐变量，写出完全数据的对数似然函数

我们设想数据是这样产生的：

首先依概率 $α_{k}$ 选择第 $k$ 个高斯分布模型 $ϕ (x | θ_{k})$ ；然后依第 $k$ 个分模型的概率分布 $ϕ (x | θ_{k})$ 生成观测数据 $x_{j}$ 。这时观测数据 $x_{j}$ ，是已知的；反映观测数据 $x_{j}$ 来自第 $k$ 个分模型的数据是未知的， $k = 1, 2, . . ., K$ ，以隐变量 $γ_{j k}$ 表示，其定义如下：

γ_{j k} = {\begin{cases} 1, 第 j 个 观 测 来 自 第 k 个 分 模 型 \\ 0, 否 则 \end{cases} j = 1, 2, . . ., N; k = 1, 2, . . ., K

有了观测数据 $x_{j}$ 及未观测数据 $γ_{j k}$ ，那么完全数据是： $(x_{j}, γ_{j 1}, γ_{j 2}, . . ., γ_{j K}) ， j = 1, 2, . . ., N$

于是，可以写出完全数据的似然函数：

p (x, γ | θ) = \prod_{j = 1}^{N} p (x_{j}, γ_{j 1}, γ_{j 2}, . . ., γ_{j K} | θ) = \prod_{k = 1}^{K} \prod_{j = 1}^{N} [α_{k} ϕ (x_{j} | θ_{k})]^{γ_{j k}} = \prod_{k = 1}^{K} α_{k}^{n_{k}} \prod_{j = 1}^{N} [ϕ (x_{j} | θ_{k})]^{γ_{j k}} = \prod_{k = 1}^{K} α_{k}^{n_{k}} \prod_{j = 1}^{N} [\frac{1}{\sqrt{2 π} σ_{k}} e x p (- \frac{(x_{j} - μ_{k})^{2}}{2 σ_{k}^{2}})]^{γ_{j k}}

式中：

n_{k} = \sum_{j = 1}^{N} γ_{j k}, \sum_{k = 1}^{K} n_{k} = N

那么，对数似然函数为：

l o g p (x, γ | θ) = \sum_{k = 1}^{K} {n_{k} l o g α_{k} + \sum_{j = 1}^{N} γ_{j k} [l o g (\frac{1}{\sqrt{2 π}}) - l o g σ_{k} - \frac{1}{2 σ_{k}^{2}} (x_{j} - μ_{k})^{2}]}

2. $E M$ 算法 $E$ 步：确定Q函数

Q (θ, θ^{(i)}) = E [l o g p (x, γ | θ) | x, θ^{(i)}] = E {\sum_{k = 1}^{K} {n_{k} l o g α_{k} + \sum_{j = 1}^{N} γ_{j k} [l o g (\frac{1}{\sqrt{2 π}}) - l o g σ_{k} - \frac{1}{2 σ_{k}^{2}} (x_{j} - μ_{k})^{2}]}} = \sum_{k = 1}^{K} {\sum_{j = 1}^{N} (E γ_{j k}) l o g α_{k} + \sum_{j = 1}^{N} (E γ_{j k}) [l o g (\frac{1}{\sqrt{2 π}}) - l o g σ_{k} - \frac{1}{2 σ_{k}^{2}} (x_{j} - μ_{k})^{2}]}

这里需要计算

E (γ_{j k} | x, θ)

，记为

\hat{γ_{j k}}

\hat{γ_{j k}} = E (γ_{j k} | x, θ) = p (γ_{j k} = 1 | x, θ) = \frac{p (γ_{j k} = 1, x_{j} | θ)}{\sum_{k = 1}^{K} p (γ_{j k} = 1, x_{j} | θ)} = \frac{p (x_{j} | γ_{j k} = 1, θ) p (γ_{j k} = 1 | θ)}{\sum_{k = 1}^{K} p (x_{j} | γ_{j k} = 1, θ) p (γ_{j k} = 1 | θ)} = \frac{α_{k} ϕ (x_{j} | θ_{k})}{\sum_{k = 1}^{K} α_{k} ϕ (x_{j} | θ_{k})}, j = 1, 2, . . ., N; k = 1, 2, . . ., K

$\hat{γ_{j k}}$ 是在当前模型参数下第 $j$ 个观测数据来自第 $k$ 个分模型的概率，称为分模型 $k$ 对观测数据 $x_{j}$ 的响应度。

将 $\hat{γ_{j k}} = E γ_{j k} 及 n_{k} = \sum_{j = 1}^{N} E γ_{j k}$ 代入之前式子即得：

Q (θ, θ^{(i)}) = \sum_{k = 1}^{K} {n_{k} l o g α_{k} + \sum_{k = 1}^{K} \hat{γ_{j k}} [l o g (\frac{1}{\sqrt{2 π}}) - l o g σ_{k} - \frac{1}{2 σ_{k}^{2}} (x_{j} - μ_{k})^{2}]}

3. $E M$ 算法 $M$ 步：

迭代的 $M$ 步是求函数 $Q (θ, θ^{(i)})$ 对 $θ$ 的极大值，即求新一轮迭代的模型参数：

θ^{(i + 1)} = a r g max_{θ} Q (θ, θ^{(i)})

求 $\hat{μ_{k}}, \hat{σ_{k}^{2}}$ 只需要将Q函数对 $μ_{k}, σ_{k}^{2}$ 求偏导并令为0，即可得到；求 $\hat{α_{k}}$ 是在 $\sum_{i = k}^{K} α_{k} = 1$ 条件下求偏导数并令为0得到的，用到拉格朗日乘子法。
结果如下：

\hat{μ_{k}} = \frac{\sum_{j = 1}^{N} \hat{γ_{j k}} x_{j}}{\sum_{j = 1}^{N} \hat{γ_{j k}}}, k = 1, 2, . . ., K \hat{σ_{k}^{2}} = \frac{\sum_{j = 1}^{N} \hat{γ_{j k}} (x_{j} - μ_{k})^{2}}{\sum_{j = 1}^{N} \hat{γ_{j k}}}, k = 1, 2, . . ., K \hat{α_{k}} = \frac{n_{k}}{N} = \frac{\sum_{j = 1}^{N} \hat{γ_{j k}}}{N}, k = 1, 2, . . ., K

重复以上计算，直到对数似然函数值不再有明显变化为止。

高斯混合模型参数估计的 $E M$ 算法

输入：观测数据 $x_{1}, x_{2}, . . ., x_{N}$ ，高斯混合模型

输出：高斯混合模型参数

（1）取参数的初始值开始迭代
（2）E步：依据当前模型参数，计算分模型 $k$ 对观测数据 $x_{i}$ 的响应度

\hat{γ_{j k}} = \frac{α_{k} ϕ (x_{j} | θ_{k})}{\sum_{k = 1}^{K} α_{k} ϕ (x_{j} | θ_{k})}, j = 1, 2, . . ., N; k = 1, 2, . . ., K

（3）M步：计算新一轮迭代的模型参数

\hat{μ_{k}} = \frac{\sum_{j = 1}^{N} \hat{γ_{j k}} x_{j}}{\sum_{j = 1}^{N} \hat{γ_{j k}}}, k = 1, 2, . . ., K \hat{σ_{k}^{2}} = \frac{\sum_{j = 1}^{N} \hat{γ_{j k}} (x_{j} - μ_{k})^{2}}{\sum_{j = 1}^{N} \hat{γ_{j k}}}, k = 1, 2, . . ., K \hat{α_{k}} = \frac{n_{k}}{N} = \frac{\sum_{j = 1}^{N} \hat{γ_{j k}}}{N}, k = 1, 2, . . ., K

（4）重复第（2）步和第（3）步，直到收敛。

至此，高斯混合模型的 $E M$ 算法我们就介绍完了，关于算法 $M$ 步中的计算，这里没有给出，求解方法会使用到带约束条件的拉格朗日乘子法
参考资料：
李航《统计学习方法》，徐亦达机器学习

【机器学习】EM算法在高斯混合模型学习中的应用

前言

单个高斯分布

高斯混合模型

高斯混合模型参数估计的EMEM算法

高斯混合模型参数估计的EMEM算法

相关推荐

高斯混合模型参数估计的 $E M$ 算法

高斯混合模型参数估计的 $E M$ 算法