EM估计（Expectation Maximization Algorithm）
这篇文章实际上是Ng在CS229上一个课程讲座。因为EM估计也很常用，特别是语义SLAM里面，这里根据一些粗浅的理解做一个总结备忘。
EM算法1是分为E（期望）和M（极大）的两步迭代优化算法，主要用于解决数据缺失情况下的概率估计问题。这个数据缺失其实不太好理解，有一个男女生身高的经典例子比较形象2，可以从ML到EM估计分别理解下。
1、极大似然估计（ML）
我们需要调查我们学校的男生和女生的身高分布。 假设你在校园里随便找了100个男生和100个女生。他们共200个人。将他们按照性别划分为两组，然后先统计抽样得到的100个男生的身高。假设他们的身高是服从正态分布的。但是这个分布的均值 \(\mu\) 和方差 \(\sigma^2\) 我们不知道，这两个参数就是我们要估计的。记作 \(\theta=[\mu,\sigma]^T\)。 这是一个典型的极大似然估计问题，我们知道100个男生，100个女生和他们的身高，我们想要估计男生的身高分布均值方差，女生的身高分布均值方差，都可以用极大似然估计来解决。
我们已知的有两个：样本服从的分布模型、随机抽取的样本；我们未知的有一个：模型的参数。根据已知条件，通过极大似然估计，求出未知参数。总的来说：极大似然估计就是用来估计模型参数的统计学方法。
因为我们知道哪个样本是男生哪个样本是女生，所以这一问题就分解成分别估计男生和女生的\(\theta=[\mu,\sigma]^T\)即可。那么以男生为例，设样本集 \(X={x_{1},x_{2},…,x_{N}}\)，其中 \(N=100\)，概率密度函数 \(p(x_{i}|\theta)\)，首先写出似然函数也就是所有男生概率的乘积：
\(l(\theta)=L\left(x_{1}, \cdots, x_{n} ; \theta\right)=\prod_{i=1}^{n} p\left(x_{i} ; \theta\right), \theta \in \Theta \tag{1}\)
通常转为求解对数似然函数：
\(H(\theta)=\ln L(\theta)=\ln \prod_{i=1}^{n} p\left(x_{i} ; \theta\right)=\sum_{i=1}^{n} \ln p\left(x_{i} ; \theta\right) \tag{2}\)
极大似然估计就是求似然函数最大值（或者对数似然函数最小值）：
\(\hat{\theta}=\arg \max _{\theta} l(\theta)=\arg \max _{\theta} \prod_{i=1}^{N} p\left(x_{i} | \theta\right) \tag{3}\)
这个极值自然就是通过求导得到：
\(\frac{d l(\theta)}{d \theta}=0 或者 \frac{d H(\theta)}{d \theta}=\frac{d \ln l(\theta)}{d \theta}=0 \tag{4}\)
特别的，对于正态分布样本 \(N\left(\mu, \sigma^{2}\right)\)，其似然函数为：
\(L\left(\mu, \sigma^{2}\right)=\prod_{i=1}^{N} \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{\left(x_{i}-\mu\right)^{2}}{2 \sigma^{2}}}=\left(2 \pi \sigma^{2}\right)^{-\frac{n}{2}} e^{-\frac{1}{2 \sigma^{2}} \sum_{l=1}^{n}\left(x_{i}-\mu\right)^{2}} \tag{5}\)
对数为：
\(\ln L\left(\mu, \sigma^{2}\right)=-\frac{n}{2} \ln (2 \pi)-\frac{n}{2} \ln \left(\sigma^{2}\right)-\frac{1}{2 \sigma^{2}} \sum_{i=1}^{n}\left(x_{i}-\mu\right)^{2} \tag{6}\)
求导得方程组：
\(\left\{\begin{array}{l}{\frac{\partial \ln L\left(\mu, \sigma^{2}\right)}{\partial \mu}=\frac{1}{\sigma^{2}} \sum_{i=1}^{n}\left(x_{i}-\mu\right)} \\ {\frac{\partial \ln L\left(\mu, \sigma^{2}\right)}{\partial \sigma^{2}}=-\frac{n}{2 \sigma^{2}}+\frac{1}{2 \sigma^{4}} \sum_{i=1}^{n}\left(x_{i}-\mu\right)^{2}=0}\end{array}\right. \tag{7}\)
则似然方程的唯一解为：
\(\left\{\begin{array}{l}{\mu^{*}=\overline{x}=\frac{1}{n} \sum_{i=1}^{n} x_{i}} \\ {\sigma^{* 2}=\frac{1}{n} \sum_{i=1}^{n}\left(x_{i}-\overline{x}\right)^{2}}\end{array}\right. \tag{8}\)
这个就是它的最大值点。也就是通过极大似然估计给出的符合当前样本的参数估计。
2、Jensen不等式
对于一个凸函数 \(f(X)\)（凸函数的定义为：对于X为实变量，\(f^{\prime \prime}(x) \geq 0 \text { (for all } x \in \mathbb{R} )\)，对于X为向量，），其中X是随机变量则有如下不等式成立：
\(E[f(X)] \geq f(E X) \tag{9}\)
其中等号成立的条件当且仅当\(X=EX\)，也就是随机变量为一个常数c：\(X=c\)。
Jensen不等式直观理解可以参见下图：
图中，实线f是凸函数，X是随机变量，有0.5的概率是a，有0.5的概率是b。（就像掷硬币一样）。X的期望值就是a和b的中值了，图中可以看到 \(E[f(X)] \geq f(E X)\) 成立。
对于f是凹函数的情况，则正相反：\(E[f(X)] \leq f(E X)\)。
3、最大期望估计（EM）
对比极大似然估计，EM估计理解起来非常简单，例如对于上述男女生的例子，我们现在有200个男女生的样本，但是我们不知道哪个样本属于男生哪个样本属于女生（信息缺失）。
3.1 EM估计公式推导
样本集 \(X={x_{1},x_{2},…,x_{m}}\)，包含 m 个独立的样本；每个样本 \(x_i\) 对应的类别 \(z_i\) 是未知的（即上文中每个样本属于哪个分布是未知的）；我们需要估计概率模型 \(p(x,z)\) 的参数 \(\theta\)，即需要找到适合的 \(\theta\) 和 \(z\) 让 \(L(\theta)\) 最大。根据上文极大似然估计中的似然函数取对数所得 \(logL(\theta)\)，可以得到如下式：
\(\ell(\theta) =\sum_{i=1}^{m} \log p(x ; \theta) \tag{10}\)
\(=\sum_{i=1}^{m} \log \sum_{z} p(x, z ; \theta) \tag{11}\)
其中第二项 \(z\) 的求和在样本缺失情况下不容易做到。但是如果 \(z\) 的概率确定后，这个求和就很容易了，EM估计就是分成两个步骤来迭代实现了这样一个估计过程。
定义 \(Q_i(z)\) 表示 \(z\) 的概率，显然 \(\sum_{z} Q_{i}(z)=1, Q_{i}(z) \geq 0\)。将上述公式展开：
\(\sum_{i} \log p\left(x^{(i)} ; \theta\right) =\sum_{i} \log \sum_{z^{(i)}} p\left(x^{(i)}, z^{(i)} ; \theta\right) \tag{12} \)
\(=\sum_{i} \log \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \tag{13} \)
\(\geq \sum_{i} \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \tag{14}\)
式 (13) 就是分子分母同乘概率 \(Q_i(z)\)，但是显然和的对数求导并不好求；式 (14) 根据 Jensen不等式而来，其中 \(f(x)=log(x)\) 是一个凹函数，同时设 \(X=\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)}\)，则 \(EX=\sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right)\left[\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)}\right]\)，根据 Jenson 不等式 \(f(EX) \geq E[f(X)]\)，则有：
\(\log \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \geq \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \tag{15}\)
再对i求和就得到了我们想要的公式。
根据之前的结论，仅当随机变量为常数时等号成立则有：
\(\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)}=c \tag{16}\)
显然 c 与 \(z_i\) 无关，则说明 \(Q_{i}\) 只和 \(p\left(x^{(i)}, z^{(i)} ; \theta\right)\) 有关，我们可以用后者作为前者的概率表达：
\(Q_{i}\left(z^{(i)}\right) \propto p\left(x^{(i)}, z^{(i)} ; \theta\right) \tag{17}\)
则显然可以推导：
\(\begin{aligned} Q_{i}\left(z^{(i)}\right) &=\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{\sum_{z} p\left(x^{(i)}, z ; \theta\right)} \\ &=\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{p\left(x^{(i)} ; \theta\right)} \\ &=p\left(z^{(i)} | x^{(i)} ; \theta\right) \tag{20} \end{aligned} \)
至此，我们推出了在固定其他参数 \(\theta\) 后，\(Q_{i}\left(z^{(i)}\right)\) 的计算公式就是后验概率，解决了 \(Q_{i}\left(z^{(i)}\right)\) 如何选择的问题。这一步就是E步，建立 \(logL(\theta)\) 的下界。接下来的M步，就是在给定 \(Q_{i}\left(z^{(i)}\right)\) 后，调整\(\theta\)，去极大化 \(logL(\theta)\) 的下界，是一个典型的极大似然估计过程（在固定 \(Q_{i}\left(z^{(i)}\right)\) 后，下界还可以调整的更大）。
如果继续用上面形象一点的例子就是：
E步：估计出男生和女生的概率分布；
M步：在上述概率分布的前提下，用极大似然估计得到当前最优的参数。
反复迭代直到收敛。
3.2 EM估计的基本步骤
基于以上推导，EM估计的典型迭代步骤如下：
循环重复直到收敛 {
(E-步骤) 对于每一个i，计算概率： \(Q_{\mathrm{i}}\left(z^{(\mathrm{i})}\right) :=p\left(z^{(i)} | x^{(\mathrm{i})} ; \theta\right)\)
(M-步骤) 使用极大似然估计最优参数：\( \theta :=\arg \max _{\theta} \sum_{i} \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \)
}
3.3 EM估计有效性的证明
我们需要证明上述的步骤一定可以得到收敛的结果，也就是 \(\ell\left(\theta^{(t)}\right) \leq \ell\left(\theta^{(t+1)}\right)\) 我们整体的极大似然估计在每次迭代后一定单调增加，那么我们一定可以在一定次数的迭代后极大似然估计不再增加也就得到了整个系统的概率最大值的参数估计，这一参数就同时包含了极大似然估计 \(\theta\) 以及缺失的参数 \(z\)。
证明如下：
对于给定 \(\theta^{(t)}\) 我们在 E 步骤选定概率函数如下：
\(Q_{i}^{(t)}\left(z^{(i)}\right) :=p\left(z^{(i)} | x^{(i)} ; \theta^{(t)}\right) \tag{21}\)
根据之前证明，这一概率函数 \(Q_{i}^{(t)}\left(z^{(i)}\right)\) 使得前述公式中的 Jensen 不等式等号成立也就是：
\(\ell\left(\theta^{(t)}\right)=\sum_{i} \sum_{z^{(i)}} Q_{i}^{(t)}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta^{(t)}\right)}{Q_{i}^{(t)}\left(z^{(i)}\right)} \tag{22}\)
则有如下推导：
\(\ell\left(\theta^{(t+1)}\right) \geq \sum_{i} \sum_{z^{(i)}} Q_{i}^{(t)}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta^{(t+1)}\right)}{Q_{i}^{(t)}\left(z^{(i)}\right)} \tag{23}\)
\(\geq \sum_{i} \sum_{z^{(i)}} Q_{i}^{(t)}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta^{(t)}\right)}{Q_{i}^{(t)}\left(z^{(i)}\right)} \tag{24}\)
\(=\ell\left(\theta^{(t)}\right) \tag{25}\)
其中式 (23) 与前述 3.1 中的公式推导一样，只是 \(\theta^{(t)}\) 换成了 \(\theta^{(t+1)}\)。式 (24) 是因为 \(\theta^{(t+1)}\) 是对公式 (22) 进行极大似然估计的结果，也就是 \(\theta^{(t+1)}\) 是使得公式中 \(\theta^{(t)}\) 最大的那个值：
\(\theta^{t+1} = \arg \max _{\theta} \sum_{i} \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \tag{26}\)
显然用 \(\theta^{(t+1)}\) 代替 \(\theta^{(t)}\) 是递增的。
根据上面的证明，可以知道 \(\ell\left(\theta^{(t)}\right) \leq \ell\left(\theta^{(t+1)}\right)\)，那么我们一定能够通过 EM估计得到最优的结果。
如果我们定义：
\(J(Q, \theta)=\sum_{i} \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \tag{27}\)
则EM估计可以看过是J的坐标上升法：，E步固定 \(\theta\)，优化 \(Q\)，M步固定 \(Q\) 优化 \(\theta\)。
参考资料
论文材料
- [pdf-embedder url="http://cdn.liuxiao.org/wp-content/uploads/2019/05/cs229-notes8.pdf"]
- [pdf-embedder url="http://cdn.liuxiao.org/wp-content/uploads/2019/05/Dempster77.pdf"]
- [pdf-embedder url="http://cdn.liuxiao.org/wp-content/uploads/2019/05/Likelihood_EM_HMM_Kalman.pdf"]
