EM算法

1. EM算法的基本思想和主要步骤

有一些优化问题、参数估计问题，往往是没有解析解的，为了得到一个近似的最优解，只能通过交替求解的方式，一步步的迭代得到近似的最优值。主要内容参考《统计学习方法》，仅作为自己学习笔记。

EM算法的主要步骤

1. 选择模型的初值，开始迭代；
2. E步：记\theta_{n}为第n次迭代参数\theta的估计值，求出Q函数的期望值Q(\theta,\theta_{n});
3. M步：求使得Q(\theta,\theta_{n})取得最大值的\theta，令\theta=\theta_{n+1}作为第n+1次迭代中求Q函数的输入参数。
4. 重复第2步和第3步，直到收敛。
   说明：这里Q函数暂且可以当作一个迭代的媒介，它是EM算法的核心。Q函数是什么东西？具体是如何推导？请看接下来的内容。

2. EM算法的导出

问题模型描述

有些参数估计问题中，往往存在显变量Y、隐变量Z以及模型的参数\theta，Y表示可观测的数据，Z表示未观测的数据，Y和Z在一起称为完全数据，Y被称为不完全数据。那么可以得出完全数据和不完全数据的似然函数、对数似然函数：
完全数据的似然函数：P_{\theta}(Y,Z)
完全数据的对数似然函数：logP_{\theta}(Y,Z)
不完全数据的似然函数：P_{\theta}(Y)
不完全数据的对数似然函数：logP_{\theta}(Y)

由最大似然估计导出EM算法

简单来想，我们是怎么来估计出模型的参数呢？那就是根据可观测变量来进行分析，为啥呢？因为不可观测变量是我们不知道的啊！那么具体怎么利用可观测变量来进行估计呢？答案就是从他的最大似然函数入手！
假设可观测变量Y=\{y_1,y_2,...,y_N\}的最大似然函数：
P_{\theta}(Y)=\prod_{j=1}^{N}P_{\theta}(y_j)
根据贝叶斯全概率公式P(Y)=\sum_z P(Y|Z)P(Z)可继续推导：
\prod_{j=1}^{N}P_{\theta}(y_j)=\prod_{j=1}^{N}\sum_{z}P_{\theta}(y_j|z)P_{\theta}(z)
似然函数是连乘的形式不好计算，因此我们定义对数似然函数：
L(\theta)=log\sum_{j=1}^{N}P_{\theta}(y_j)=log\sum_{j=1}^{N}\sum_{z}P_{\theta}(y_j|z)P_{\theta}(z)
为了方便推导，简写为：
L(\theta)=log\sum_{z}P_{\theta}(Y|Z)P_{\theta}(Z)

有了对数似然函数以后，就可以按照套路来估计参数\theta。什么意思呢？意思就是我们通过似然函数计算出一个\theta，使得L(\theta)取得最大值。即：
\hat{\theta}=arg\max_{\theta}L(\theta)
对于这个问题是求不出解析解的，所以需要迭代求解的方法。而EM算法就是用于求解这种问题的一个迭代求解方法。

EM算法的推导

EM算法怎么就可以迭代求解得出最优值呢？其关键之处就在于两个重要的步骤：(1)E步，求期望;(2)M步，最大化。

E步：

有一个迭代的前提条件，如果最大似然函数L(\theta_{n}) < L(\theta_{n+1})，那么就可以继续迭代求出更优的\theta。否则就找到了局部或全局的最优解。
那么满足前提条件L(\theta_{n}) < L(\theta_{n+1})的时候怎么进行E步的求期望呢？求的是什么期望呢？怎么就需要求这个期望呢？
这些问题乍一看很让人头大，实际上还需要从数学的角度去考虑，“期望”是什么？概率论里的期望不就\sum_{i=1}^{N}x_i \cdot P(x_i)吗？
那么EM算法里的期望是什么呢？它指的是一个下边界函数Q的期望。Q函数怎么来的呢？请看推导：
我们为了得到增量，需要计算：
L(\theta)-L(\theta_n)=\sum_{z}P_{\theta}(Y|Z)P_{\theta}(Z)-\sum_{z}P_{\theta_n}(Y|Z)P_{\theta_n}(Z)
上面的式子中，\theta是变量，\theta_n是常量，即第n次迭代中的\theta值。因此我们可以根据贝叶斯全概率公式将其简化的得到：
L(\theta)-L(\theta_n)=log\sum_{z}P_{\theta}(Y|Z)P_{\theta}(Z)-log\sum_{z}P_{\theta_n}(Y)
对于上式中的两项，我们需要通过一些数学技巧来凑到一起，当然这些奇怪的想法是那些牛批的大佬想出来的，打死我也想不出这么巧妙的推导。
第一项：
log\sum_{z}P_{\theta}(Y|Z)P_{\theta}(Z)=log\sum_{z}P_{\theta_n}(Y|Z)\cdot \frac{P_{\theta}(Y|Z)P_{\theta}(Z)}{P_{\theta_n}(Y|Z)}
注意到这是一个关于\theta的函数，因此我们根据凸函数的詹森不等式：f(EX) \geqslant E[f(x)]可以把函数log丢进去，继续推导：
\begin{aligned} log\sum_{z}P_{\theta}(Y|Z)P_{\theta}(Z) & =log\sum_{z}P_{\theta_n}(Y|Z)\cdot \frac{P_{\theta}(Y|Z)P_{\theta}(Z)}{P_{\theta_n}(Y|Z)}\\ & \geqslant \sum_{z}P_{\theta_n}(Y|Z)\cdot log\frac{P_{\theta}(Y|Z)P_{\theta}(Z)}{P_{\theta_n}(Y|Z)}\\ \end{aligned}
这样我们就得到了第一项，第二项还是第二项，没有变，把上面的式子整理下来：
\begin{aligned} L(\theta)-L(\theta_n)&=\sum_{z}P_{\theta}(Y|Z)P_{\theta}(Z)-\sum_{z}P_{\theta_n}(Y|Z)P_{\theta_n}(Z)\\ &= log\sum_{z}P_{\theta}(Y|Z)P_{\theta}(Z) -log\sum_{z}P_{\theta_n}(Y)\\ &=log\sum_{z}P_{\theta_n}(Y|Z)\cdot \frac{P_{\theta}(Y|Z)P_{\theta}(Z)}{P_{\theta_n}(Y|Z)}-log\sum_{z}P_{\theta_n}(Y)\\ &\geqslant \sum_{z}P_{\theta_n}(Y|Z)\cdot log\frac{P_{\theta}(Y|Z)P_{\theta}(Z)}{P_{\theta_n}(Y|Z)}-log\sum_{z}P_{\theta_n}(Y)\\ &=\sum_{z}P_{\theta_n}(Y|Z)\cdot log\frac{P_{\theta}(Y|Z)P_{\theta}(Z)}{P_{\theta_n}(Y|Z)P_{\theta_n}(Y)} \end{aligned}
把L(\theta_n)移到不等号右边，令
B(\theta,\theta_n)=L(\theta_n)+\sum_{z}P_{\theta_n}(Y|Z)\cdot log\frac{P_{\theta}(Y|Z)P_{\theta}(Z)}{P_{\theta_n}(Y|Z)P_{\theta_n}(Y)}
得到了下边界函数B(\theta,\theta_n)，它具有如下两个性质：
性质1. L(\theta_n)=B(\theta_n,\theta_n)
性质2. L(\theta) \geqslant B(\theta,\theta_n)
性质1说明当\theta=\theta_n时，L(\theta)与下边界函数B(\theta,\theta_n)有一个交点；
性质2说明B(\theta,\theta_n)总是小于L(\theta)。
因此，根据性质2，当我们找到一个\theta，能使B(\theta,\theta_n)增大时，L(\theta)也会增大。为了使得L(\theta)尽可能增大，我们就要最大化B(\theta,\theta_n)，并找出能使其最大化的\theta_{n+1}。这就引出了M步。
TIPS:上面为什么要让L(\theta)尽可能增大呢？原因是当L(\theta)取得极大值时，L(\theta)-L(\theta_n)=0时，迭代就可以结束了，就可以得到最优值了。

M步：

上面得出M步需要最大化B(\theta,\theta_n)，并得出\theta_{n+1}，即：
\hat{\theta}_{n+1}=arg\max_{\theta}B(\theta,\theta_n)

导出Q函数

实际上，对于arg\max_{\theta}B(\theta,\theta_n)可以作一下简化，去掉与\theta无关的常数项：
\begin{aligned} \hat{\theta}_{n+1}&=arg\max_{\theta}B(\theta,\theta_n)\\ &=arg\max_{\theta}\big[L(\theta_n)+\sum_{z}P_{\theta_n}(Y|Z)\cdot log\frac{P_{\theta}(Y|Z)P_{\theta}(Z)}{P_{\theta_n}(Y|Z)P_{\theta_n}(Y)}\big]\\ &=arg\max_{\theta}\big[L(\theta_n)+\sum_{z}P_{\theta_n}(Y|Z)\cdot logP_{\theta}(Y|Z)P_{\theta}(Z)-\sum_{z}P_{\theta_n}(Y|Z)\cdot P_{\theta_n}(Y|Z)P_{\theta_n}(Y)\big]\\ &=arg\max_{\theta}\sum_{z}P_{\theta_n}(Y|Z)\cdot logP_{\theta}(Y|Z)P_{\theta}(Z)\\ &=arg\max_{\theta}Q(\theta,\theta_n)\\ \end{aligned}
这个Q函数实际上就是在最大化过程，简化下边界函数B(\theta,\theta_n)得到的结果。通过上式的最大化，我们就导出了\theta_{n+1}，这相当于一次迭代。
用\theta_{n+1}可以求Q(\theta,\theta_{n+1})，并通过Q(\theta,\theta_{n+1})导出\theta_{n+2}；
用\theta_{n+2}求Q(\theta,\theta_{n+2})，并通过Q(\theta,\theta_{n+2})导出\theta_{n+3}；
...........
最后达到收敛条件，停止迭代，此时求到的\theta^*就是我们最终用EM算法估计出来的模型参数。
完毕！