彻底理解机器学习

阿南带你理解EM算法

以为隐藏起来我就找不到你了？

Posted by Nam on February 27, 2020

前言

前面讲了那么多的机器学习方法，其都是基于模型确定的情况下通过训练得到决策函数。想必各位看官也有些厌倦了~

今天给大家分享的就是我学了这么久，觉得最不好理解的一个算法 —— EM算法。

EM 算法

老规矩，先解释一下 EM 算法 这个名字是咋来的。 EM算法 中的 $E$ 是 Expectation 期望的缩写，$M$ 是 Maximization 极大化的缩写。所以 EM 算法 每次计算分为两步进行：$E$ 步求期望，$M$ 步求极大。

那么问题来了，EM 算法 有什么独到之处呢？为啥子要学这个算法？

柴犬

既然你诚心诚意地发问了，那我就非常认真地告诉你：在现实生活中，事物的组成往往是多元化的，比方说在学校里面随机抽选100个学生，这100个学生中也是有男有女而不会说是只有一种性别（女子学校除外哈）。

试想这么一个问题，我们想知道抽出的这一百个学生中男生和女生的平均身高，借此作为全校男生和女生身高的无偏估计。但是问题来了，当时记录身高的人是个马大哈，他只记录了这一百个学生的身高却没有记录对应的性别。那再这种情况下我们应该怎么做呢？

重新抽100个学生再次进行记录吗？ NO NO NO NO！ 这个不符合我们高贵的气质和深厚的技术底蕴。阿南告诉你，我们只要用今天的主角 EM 算法 就可以解决这种已知样本分布的种类和数量，但是不知道每个样本 $x$ 来自何种分布 $dist$，求解各个分布参数的问题。

OK! 现在假设男女生身高服从正态分布。（万物皆可正态！）

\[\phi (x|\theta_k) = \frac{1}{\sqrt{2\pi}\sigma_k} \times e^{ - \frac{(x-\mu_k)^2}{2\sigma^2_k}} \tag{1}\]

其中 $\theta_k$ 是待求的参数 $\theta_k=(\mu_k,\sigma^2_k)$，密度函数图像如下所示：

EM

（图画滴有点丑，大家将就看哈~）

接下来我们的目标就是分别求出上图中两种分布的均值分别作为男生和女生平均身高的 无偏估计。

问题是，我们这100个身高数据中不知道每个数据对应的性别，咋么求对应分布的均值捏？这时候 EM 算法 站出来说：”我可以，我能行！“

我可以

EM 算法 中将不可观测的变量 $Z$ 称之为 隐变量。

用 $X$ 表示随机变量的数据，$Z$ 表示隐随机变量的数据。在上述身高问题中，身高样本 $x_i$ 来自何种分布 $dist$ 就是不可观测的隐变量 $Z$。$X$ 和 $Z$ 连在一起称为完全数据，单独观测数据称为不完全数据。假设给定观测数据 $X$ ，其概率分布是 $P(X\mid \theta)$，其中 $\theta$ 是我们需要估计的模型参数。

那么不完全数据 $X$ 的似然函数是 $P(X\mid \theta)$，对数似然函数 $L(\theta)=\log P(X\mid \theta)$。假设 $X$ 和 $Z$ 的联合概率分布是 $P(X,Z\mid \theta)$，那么完全数据的似然函数是 $\log P(X,Z\mid \theta)$

在上述身高问题中，由于数据来 $x$ 自两个参数不同的高斯分布，我们称之为 混合高斯分布 (Gaussian Mixture Model, GMM)，具有如下形式的概率分布模型：

\[P(x|\theta) = \sum_{k=1}^K \pi_k \phi(x|\theta_k) \tag{2}\]

其中，$\pi_k$ 是每个分布出现的概率，$\pi_k \ge 0,\ \sum_{k=1}^K \pi_k = 1;\ \phi(x\mid \theta_k)$ 是高斯分布的密度函数（见式 $1$)，$\theta_k = (\mu_k,\sigma_k^2)$。

这里给出 EM算法 的核心思想：EM算法通过迭代求解观测数据的对数似然函数 $L(\theta)=\log P(X\mid \theta)$ 的极大化，实现对待求参数 $\theta$ 的极大似然估计。 即极大化：

\[\begin{aligned} L(\theta) & = \log P(X|\theta) \\[2ex] & = \log \sum_{Z} P(X,Z|\theta) \\ & = \log [\sum_{Z} P(X|Z,\theta)P(Z|\theta)] \end{aligned} \tag{3}\]

仔细观察可以发现，极大化式 $(3)$ 的难点在于：由于未观测数据的存在导致对数中出现加和 $\sum$ 导致优化求解的过程非常困难。事实上，EM算法 是通过迭代逐步近似极大化 $L(\theta)$ 的。假设在第 $i$ 次迭代以后 $\theta$ 的估计值是 $\theta^{(i)}$。我们希望新估计值 $\theta$ 能够使 $L(\theta)$ 增加，即 $L(\theta)>L(\theta^{(i)})$，并逐步达到极大值。

在这里我们利用 $Jensen$ 不等式 可以得到 $L(\theta)$ 式的下界：

\[\begin{aligned} L(\theta) & = \log (\sum_{Z} P(X|Z,\theta)P(Z|\theta)) \\ & \ge \sum_{Z} \{P(Z|X,\theta^{(i)}) \log \frac{P(X|Z,\theta)P(Z|\theta)}{P(Z|X,\theta^{(i)})}\} \\ \end{aligned} \tag{4}\]

注：因为 $log$ 函数是单调递增凸函数，根据 $Jensen$不等式，期望的 $log$ 大于 $log$ 的期望。形象的解释如下图：

Jensen

令

\[B(\theta, \theta^{(i)}) = \sum_{Z} \{P(Z|X,\theta^{(i)}) \log \frac{P(X|Z,\theta)P(Z|\theta)}{P(Z|X,\theta^{(i)})}\} \tag{5}\]

则有

\[L(\theta) \ge B(\theta, \theta^{(i)}) \tag{6}\]

即函数 $B(\theta,\theta^{(i)})$ 是 $L(\theta)$ 的一个下界。因此，任何可以使 $B(\theta, \theta^{(i)})$ 增大的参数 $\theta$ 都可以使 $L(\theta)$ 增大。为了使 $L(\theta)$ 在一次步长内有尽可能大的增长，我们因选择 $\theta^{(i+1)}$ 使得 $B(\theta, \theta^{(i)})$ 达到极大，即：

\[\theta^{(i+1)} = \arg \max_{\theta} B(\theta,\theta^{(i)}) \tag{7}\]

现在我们求解 $\theta^{(i+1)}$ 的表达式。我们摒弃对极大化 $\theta$ 而言是常数的项，由式 $(5),(7)$ 我们可得：

\[\begin{aligned} \theta^{(i+1)} & = \arg \max_{\theta} (\sum_{Z} \{P(Z|X,\theta^{(i)}) \log \frac{P(X|Z,\theta)P(Z|\theta)}{P(Z|X,\theta^{(i)})}\}) \\ & = \arg \max_{\theta} \{\sum_{Z} [P(Z|X,\theta^{(i)}) \log P(X|Z,\theta)P(Z|\theta)]\} \\ & = \arg \max_{\theta} \{\sum_{Z} [P(Z|X,\theta^{(i)}) \log P(X,Z|\theta)]\} \\ & = \arg \max_{\theta} Q(\theta, \theta^{(i)}) \end{aligned} \tag{8}\]

式 $(8)$ 中我们将完全数据的对数似然函数 $\log P(X,Z\mid \theta)$ 关于在给定观测数据 $X$ 和当前参数 $\theta^{(i)}$ 下对未观测数据 $Z$ 的条件概率分布 $P(Z\mid X,\theta^{(i)})$ 的条件期望称之为 $Q$函数：

\[\begin{aligned} Q(\theta, \theta^{(i)}) & = E_Z[\log P(X,Z|\theta)|X,\theta^{(i)}] \\[2ex] & = \sum_{Z} [P(Z|X,\theta^{(i)}) \log P(X,Z|\theta)] \end{aligned} \tag{9}\]

因此我们只需找到能使 $Q$函数 极大化的参数 $\theta^{(i)}$ 即完成一次迭代操作。 上述操作再次表明了 EM算法 是通过不断求解下界的极大化逼近求解对数似然函数的极大化算法。

EM1

（图片来自李航老师的《统计学习方法》）

OK！这么繁杂滴概念介绍完了，现在我们就用 $EM$算法 求解 高斯混合模型(GMM):

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

    我们假设我们的数据集 $x_i,i=1,2,…,N$，是这样产生的：首先根据概率 $\pi_k$ 选择第 $k$ 个高斯分布模型 $\phi(x\mid \theta_k)$。然后按照第 $k$ 个分模型的概率分布 $\phi(x\mid \theta_k)$ 生成观测数据 $x_i$。

    但是现在问题就是我们不知道观测数据 $x_i$ 来自哪个分模型。因此在此处我们定义：观测数据 $x_i$ 来自第 $k$ 个分模型为隐变量 $Z$：

    \[Z_{ik} = \begin{cases} 1, \qquad 第 i 个观测来自第 k 个分模型 \qquad i = 1,2,...,N;\\[2ex] 0, \qquad Otherwise \qquad \qquad \qquad \ \ \ \ \qquad k=1,2,...,K \\ \end{cases} \tag{10}\]

    $Z_{ij}$ 是 $0-1$ 随机变量。

    有了观测数据 $x_i$ 及未观测数据 $Z_{ij}$，那么完全数据就是：

    \[(x_i,Z_{i1},Z_{i2},...,Z_{iK}), i=1,2,...,N\]

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

    \[\begin{aligned} P(X,Z|\theta) & = \prod_{i=1}^N P(x_i,Z_{i1},Z_{i2},...,Z_{iK}|\theta) \\ & = \prod_{k=1}^K \prod_{i=1}^N [\pi_k \phi(x_i|\theta_k)]^{Z_{ik}} \\ & = \prod_{k=1}^K \pi_k^{n_k} \prod_{i=1}^N [\phi(x_i|\theta_k)]^{Z_{ik}} \\ & = \prod_{k=1}^K \pi_k^{n_k} \prod_{i=1}^N [\frac{1}{\sqrt{2\pi}\sigma_k} \times e^{ - \frac{(x-\mu_k)^2}{2\sigma^2_k}}]^{Z_{ik}} \end{aligned} \tag{11}\]

    式 $(11)$ 中，$n_k = \sum_{i=1}^N Z_{ik},\ \ \sum_{k=1}^K n_k = N$。

    取对数，得到完全数据的对数似然函数：

    \[\log P(X,Z|\theta) = \sum_{k=1}^K \{n_k \log \pi_k + \sum_{i=1}^N Z_{ik} [\log (\frac{1}{\sqrt{2\pi}})-\log \sigma_k - \frac{1}{2\sigma^2_k}(x_i-\mu_k)^2] \} \tag{12}\]

  2. EM算法的 $E$ 步：确定 $Q$ 函数 \(\begin{aligned} Q(\theta,\theta^{(i)}) & = E[\log P(X,Z|\theta)|X,\theta^{(i)}] \\[2ex] & = E\{\sum_{k=1}^K\{n_k\log \pi_k + \sum_{j=1}^N Z_{ik} [\log (\frac{1}{\sqrt{2\pi}})- \log \sigma_k - \frac{1}{2\sigma^2_k}(x_i-\mu_k)] \} \} \\[2ex] & = \sum_{k=1}^K \{\sum_{i=1}^N \{E(Z_{ik})\log \pi_k\} +\sum_{i=1}^N (E_{Z_{ik}}) [\log (\frac{1}{\sqrt{2\pi}})- \log \sigma_k - \frac{1}{2\sigma^2_k}(x_i-\mu_k)] \} \end{aligned} \tag{13}\)

    观察上式发现存在除了待求解参数 $\theta$ 以外还存在未知变量 $E(Z_{ik})$。

    因此在 $E$ 步中我们需要根据当前拥有的参数 $\theta^{(i)}$ 先计算隐变量 $Z_{ik}$ 的数学期望 $E(Z_{ik}\mid X,\hat \theta)$，记为 $\hat Z_{ik}$

    \[\begin{aligned} \hat Z_{ik} & = E(Z_{ik}\mid X,\theta) = P(Z_{ik}=1|X,\theta) \\[2ex] & = \frac{P(Z_{ik}=1,x_i|\theta)}{\sum_{k=1}^K P(Z_{ik}=1,x_i|\theta)} \\[2ex] & = \frac{P(x_i | Z_{ik}=1,\theta)P(Z_{ik}=1|\theta)}{\sum_{k=1}^K P(x_i | Z_{ik}=1,\theta)P(Z_{ik}=1|\theta)} \\[2ex] & = \frac{\pi_k \phi(x_i|\theta_k)}{\sum_{k=1}^K \pi_k \phi(x_i|\theta_k)},\qquad i = 1,2,...,N; \qquad k=1,2,...,K \end{aligned} \tag{14}\]

    $\hat Z_{ik}$ 是在当前模型参数下第 $i$ 个观测数据来自第 $k$ 个分模型的概率，称之为分模型 $k$ 对观测数据 $x_i$ 的响应度。

    小伙伴们可能会有疑问，之前的隐变量 $Z_{ik}$ 是一个指示函数 $sign$，为啥到这里求了个期望之后就变成了概率了咧？答案是这样的：之前的确是非 $0$ 即 $1$ 的指示函数，但是其期望 $\hat Z_{ik}$ 根据期望的定义就是：

    \[\begin{aligned} E(Z_{ik}) = \hat Z_{ik} & = 1\times P(Z_{ik} =1) + 0\times P(Z_{ik} =0) \\[2ex] & = P(Z_{ik} =1) \end{aligned} \tag{15}\]

    经过式 $(15)$ 的推导，$\hat Z_{ik}$ 显然是在当前模型参数下第 $i$ 个观测数据来自第 $k$ 个分模型的概率。

    将 $\hat Z_{ik}=E(Z_{ik})$ 及 $n_k=\sum_{j=1}^N E(Z_{ik})$ 带入 $Q$函数 的表达式 $(13)$，即得：

    \[Q(\theta,\theta^{(i)}) = \sum_{k=1}^K \{ n_k \log \pi_k +\sum_{i=1}^N \hat Z_{ik} [\log (\frac{1}{\sqrt{2\pi}})- \log \sigma_k - \frac{1}{2\sigma^2_k}(x_i-\mu_k)] \} \tag{16}\]

    此时 $Q(\theta,\hat \theta)$ 就只剩下 $\theta$ 一个变量了，另外的 $\hat \theta$ 式一个确定的值，接下来便阔以计算参数 $\theta$ 了！

  3. EM算法 的 $M$ 步：极大化 $Q(\theta,\hat \theta)$ 求解参数 $\theta$

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

    \[\theta^{(i+1)} = \arg \max_\theta Q(\theta,\theta^{(i)}) \tag{17}\]

    用 $\hat \mu,\hat \sigma^2$ 及 $\hat \pi_k$ 来表示 $\theta^{(i+1)}$ 的各参数。求解 $\hat \mu^2,\hat \sigma^2$ 只需将式 $(16)$ 分别对 $\mu,\sigma^2$ 求偏导并令其为零即可得到。求解 $\hat \pi_k$ 只需根据定义：用所有样本 $x_i$ 在第 $k$ 个分模型中出现期望的累加和 除以 所有样本在所有分模型中出现的期望——即样本总数量$N$，即可。

    \[\begin{aligned} \hat \mu_k & = \frac{\sum_{i=1}^NZ_{ik}x_i}{\sum_{i=1}^NZ_{ik}} \\[2ex] \hat \sigma^2_k & = \frac{\sum_{i=1}^NZ_{ik}(x_i-\mu_k)^2}{\sum_{i=1}^NZ_{ik}} \\[2ex] \hat \pi_k & = \frac{n_k}{N} = \frac{\sum_{i=1}^NZ_{ik}}{N} \end{aligned} \tag{18}\]

    重复以上的计算过程，直至对数似然函数 $(12)$ 的值不再有明显变化时我们就得到所有 $K$ 个模型对应的参数值 $\theta_k$ 以及每个分模型出现的概率 $\pi_k$。


总结

  1. EM算法 的结果与初值的选择有关，不同的初值可能得到不一样的参数估计值。
  2. EM算法 只能得到局部最优解而不能得到全局最优解。
  3. EM算法 时通过迭代极大化观测数据的对数似然函数 $L(\theta)=\log P(X\mid \theta)$ 去实现极大似然估计。
  4. EM算法 的构建最重要的是定义 $Q$函数。因为 EM算法 极大化似然函数的方法就是在每次的迭代过程中通过极大化 $Q$函数 来间接地增大对数似然函数。

好滴， 今天的 EM算法 小课堂就到这里啦！

谢谢大家 :)

CATALOG
    FEATURED TAGS
    彻底理解机器学习