1. 最大似然估计到EM算法
之前在介绍《最大似然估计》中提到，MLE估计的目标是产生样本集合可能性最大的参数。其核心思想是在已知结果情况下，然后寻求使该结果出现的可能性最大的条件，并以此作为估计值。假设：要根据iris数据中的花蕊长度来估计花瓣的长度，我们可以利用最小二乘来估计参数。根据极大似然估计的定义，我们要获得产生样本可能性最大的参数\(\theta\)，则似然函数可以记为：
\[\ell(\theta)=\prod^m_{i=1} \frac{1}{\sqrt{2\pi}\sigma} exp(- \frac{(y^{(i)}-\theta x)^2}{2 \sigma})\]由于线性回归服从独立同分布的假设，因此，我们可以忽略\(\sigma\)的影响，直接估计\(\theta\)。
但是，我们假设如果知道不仅花蕊的长度与花瓣的长度有关，同时，花的种类对花瓣的长度有影响，即不同种类花的花瓣长度可能服从不同分布，我们仍然在独立同分布的假设下采用极大似然估计求解是不合适的。因为，不同种类花的花瓣与花蕊长度的误差项不服从同一正态分布，即不同分布的\(sigma\)不同，因此，我们无法忽略该参数对目标函数的影响，需要同时估计两个参数\(\theta\)和\(\sigma\)。但是假设我们如果无法获知各个样本隶属的花的类别，即无法通过参数\(\sigma\)的先验知识估计参数\(\sigma\)，那么我们如何来估计花瓣的长度呢？
总结上面的问题，我们发现，（1）我们不知道每个样本隶属于的花的类别；（2）我们也不知道花瓣与花蕊的线性关系是什么。从数学角度考虑这些问题，可以映射成：（1）抽样样本不知道是从哪个分布抽取的；（2）每个分布下的参数是什么。
可以确认的是，只有当我们知道每个样本属于哪个分布，我们才会准确的估计每个分布下花瓣与花蕊长度的线性关系。但是，各个类别之间的样本混在一起，我们无法获知哪个样本属于哪个类别，所以就无法准确估计各分布的参数。反之，只有当我们获得各个分布参数的准确估计下，我们才能计算出各样本隶属于哪个类别。因此，无法直接采用极大似然估计求解。这个问题在数学上是一个典型的循环依赖问题，即“鸡生蛋，蛋生鸡”的问题。想要解决循环依赖问题，我们可以先给定任意一个参数\(\theta\)初始值，这样，通过初始值来估计参数\(\sigma\)，然后再用估计出的参数\(\sigma\)再次估计参数\(\theta\)，通过多次的循环迭代，最终收敛到一个最优解。
EM算法（Expectation Maximization）的核心思想就是这样，假设我们要估计初始状态未知的两个参数，同时两个参数间存在依赖关系，那么我们就可以先赋予其中一个参数初始值，得到另一个参数的估计值，然后在另一个参数估计值的基础上，重新估计第一个参数，直至收敛。在上面的问题中，如果我们先给出每个花的类别下花瓣和花蕊长度间的线性关系（\(\theta\)），那么我们可以根据\(\theta\)值计算每个样本可能属于的花的类别（花瓣长度和花蕊长度的线性误差小的类别），这个过程在EM算法中属于Expectation步骤。通过Expectation步，每个样本都有所属类别，这样我们就可以将样本分成几部分，然后采用极大似然估计的思想，对每个部分中的样本重新估计参数，即估计各类别下花瓣和花蕊长度间的线性关系，这个就是EM算法中的Maximization步骤。如此迭代，直至各个参数不在发生变化即可停止迭代。
一般情况下，我们把每个样本隶属于哪种花对应的参数称为隐含变量，因为，如果没有隐含变量的话，我们可以很直观的求解待估计参数，当存在隐含变量后，本来很简单求解的问题变得复杂。因此，EM算法存在的意义就在于此，可用于包含隐含变量的参数估计。但是，按照最大似然的思想，我们每次估计会使得目标函数上升，最终达到极值，而EM算法怎样保证算法收敛？怎样保证可以获得目标函数的极值？这些问题就必须从数学角度进行推导了～
2. EM算法原理
假设给定训练集\(X={x^{(1)}, x^{(2)}, \ldots, x^{(m)}}\)，我们希望求解包含隐含变量\(Z＝{z^{(1)}, z^{(2)}, \ldots, z^{(m)}}\)的目标函数\(f(x,z;\theta)\)中的参数\(\theta\)。根据最大似然估计，我们可以获得：
\[\theta=arg max_{\theta} \ell (\theta;X)= arg max_{\theta} log \sum_Z f(X,Z;\theta)＝arg max_{\theta} \sum_{i=1}^m log \sum_{z^{(i)}} f(x^{(i)},z^{(i)};\theta)\]对于上述公式，想要求解\(log \sum_Z f(X,Z;\theta)\)形式的极值往往非常困难，如果我们可以观察到隐藏变量\(z\)的话，则可以消除\(\sum\)符号，这样求解极值就容易的多（采用梯度下降算法）。但实际上我们没法获得\(z\)的分布或者\(z\)值，只可能通过给定\(\theta\)的情况下，通过求解后验概率\(P(z｜x;\theta)\)获得参数\(z\)的估计值。
我们令\(L(\theta)=log \sum_Z P(X,Z;\theta)\)，其中\(Z\)为隐含变量。很明显，求解该log嵌套\(\sum\)形式的极值很麻烦，因此不妨引入概率分布\(Q(Z)\)，其中\(\sum_{Z} Q(Z)=1\)，且\(Q(Z) \geq 0\)，得到：
\[L(\theta)=log \sum_{Z} P(X,Z;\theta)=log \sum_{Z} Q(Z) \frac{P(X,Z;\theta)}{Q(Z)}\]a. 离散变量的期望
对于\(\frac{P(X,Z;\theta)}{Q(Z)}\)而言，隐含变量\(Z\)为离散变量，如果我们假定\(Q(Z)\)正好是离散随机变量\(Z\)的概率分布函数，且\(Y=g(Z)=\frac{P(X,Z;\theta)}{Q(Z)}\)为连续函数，那么：
\[E(Y)=E[g(Z)]=\sum_Z g(Z)Q(Z)\]如果\(Z\)是连续随机变量，那么：
\[E(Y)=E[g(Z)]=\int^{\infty}_{-\infty} g(Z)Q(Z)dZ\]所以：
\[L(\theta)=log \sum_{Z} Q(Z) \frac{P(X,Z;\theta)}{Q(Z)}=log E_Z[g(Z)]=log E_Z [\frac{P(X,Z;\theta)}{Q(Z)}]\]b. Jensen不等式
如果任意函数\(f\)为凸函数，\(X\)为随机变量，则\(E[f(X)] \geq f(E[X])\)，若上式取等号当且仅当\(X\)为常数，若任意函数\(f\)为凹函数，则不等式反向。因为log函数为严格凹函数，所以对于任意分布\(Q\)，
\[L(\theta)=log E_Z [\frac{P(X,Z;\theta)}{Q(Z)}] \geq E_Z[log \frac{P(X,Z;\theta)}{Q(Z)}]\]c. 凹函数的下界函数
对于任意凹函数\(f\)，假设我们想求\(arg max_{\theta} f(\theta)\)，对于一般优化算法而言，我们需要保证每次迭代得到的\(f(\theta_{t+1})\)一定比\(f(\theta_t)\)要大，这样才能保证优化过程呈单调逼近于极大值。
因此，如果对于任意\(t\)次迭代，我们能找到\(f(\theta)\)的下界函数\(g_{\theta_t}(\theta)\)，满足\(f(\theta) \geq g_{\theta_t}(\theta)\)，且\(f(\theta_t) \geq g_{\theta_t}(\theta_t)\)那么令\(\theta_{t+1}:=arg max_{\theta} g_{\theta_t}(\theta)\)，此时，\(g_{\theta_t}(\theta_{t+1}) ＝ max_{\theta} g_{\theta_t}(\theta_t)\),即会满足每次迭代成单调趋势：
\[f(\theta_{t+1}) \geq g_{\theta_t}(\theta_{t+1}) \geq g_{\theta_t}(\theta_t)=f(\theta_t)\]对于\(t＋1\)次迭代，我们仍然可以找到函数\(g_{\theta_{t+1}}(\theta)\)，满足\(f(\theta_{t+2}) \geq g_{\theta_{t+1}}(\theta_{t+2}) \geq g_{\theta_{t+1}}(\theta_{t+1})=f(\theta_{t+1})\)，据此计算出\(\theta_{t+2}\)。
再回到EM算法中，如果我们可以找到\(L(\theta)\)的下界函数，并且让不等式在\(\theta_t\)处取等号，即可使得每次迭代的目标函数更优。我们已知\(L(\theta) \geq E_Z[log \frac{P(X,Z;\theta)}{Q(Z)}]\)，所以可以看出\(E_Z[log \frac{P(X,Z;\theta)}{Q(Z)}]\)为下界函数，但是我们还要让下界函数在\(\theta=\theta_t\)处取等号，即要求\(\frac{P(X,Z;\theta)}{Q(Z)}\)为常数。若要使得\(Q(Z) \propto P(X,Z;\theta_t)\)，且\(\sum Q(Z) = 1\)，则根据全概率公式可知，仅当\(Q(Z)=\frac{P(X,Z;\theta_t)}{\sum_Z P(X,Z;\theta_t)}=P(Z｜X;\theta_t)\)时，满足\(Q(Z)\)与\(P(X,Z;\theta_t)\)呈正比关系。
因此，我们令\(Q(Z)=P(Z｜X;\theta_t)\)，此时：
\[L(\theta) \geq E_Z[log \frac{P(X,Z;\theta)}{Q(Z)}] ＝ E_{Z \sim Q = Z|X;\theta_t} [log \frac{P(X,Z;\theta)}{P(Z|X;\theta_t)}]\]显然，上式在\(\theta = \theta_t\)时，\(L(\theta_t) ＝ E_{Z｜X;\theta_t} [log \frac{P(X,Z;\theta_t)}{P(Z｜X;\theta_t)}] = E_{Z｜X;\theta_t} [log \sum_Z P(X,Z;\theta_t)]\)可以取等号。
综上所述，\(E_{Z｜X;\theta_t} [log \frac{P(X,Z;\theta)}{P(Z｜X;\theta_t)}]\)为\(L(\theta)\)的下界函数，我们可以通过每次迭代计算找到满足：
\[\theta_{t+1} := arg max_{\theta} E_{Z｜X;\theta_t} [log \frac{P(X,Z;\theta)}{P(Z｜X;\theta_t)}]\]d. \(\theta\)更新公式推导
找到目标函数的下界函数\(g_{\theta_t}(\theta)\)后，我们可以通过更新\(\theta_{t+1}:=arg max_{\theta} g_{\theta_t}(\theta)\)实现每次的迭代，并可以保证迭代过程呈单调变化趋势，最终收敛。
我们可以再对更新公式进行变换，使得计算更为简便：
\(\theta_{t+1} = arg max g_{\theta_t} (\theta) \\ = arg max E_{Z｜X;\theta_t} \, [log \frac{P(X,Z;\theta)}{P(Z｜X;\theta_t)}] \\ = arg max \sum_Z P(Z｜X;\theta_t) log \frac{P(X,Z;\theta)}{P(Z｜X;\theta_t)} \\ = arg max \sum_Z [P(Z｜X;\theta_t)logP(X,Z;\theta) - P(Z｜X;\theta_t)log P(Z｜X;\theta_t)] \, (\because 与 \theta无关) \\ = arg max \sum_Z [P(Z｜X;\theta_t)logP(X,Z;\theta) \\ = arg max E_{Z｜X;\theta_t} log P(X,Z;\theta)\)。
3. EM算法过程
根据上面的推导过程，我们可以定义EM算法的Expectation步和Maximization步。
Repeat until converge｛
- E－step：计算\(P(Z｜X;\theta_t)\)得到下界函数
- M－step：更新\(\theta\)
｝
4. 总结
总的来说，Expectation Maximization算法是一个比较难懂的算法，中间的推导过程较为复杂，可能学习一遍两遍很难把这个算法吃透，需要比较深入的统计学知识才能搞明白，这篇博文也写了三天才把算法原理过程写完。
EM算法是用于对含有隐含变量的目标函数做最大似然估计的一种方法。EM算法的应用特别广泛，例如，Kmeans聚类，高斯混合模型，隐马尔可夫模型等。后面的博文我会专门撰写一篇EM算法实例，用于具体说明EM算法求解过程。