机器学习基础EM算法

前言
EM算法是机器学习十大算法之一，它很简单，但是也同样很有深度，简单是因为它就分两步求解问题，

E步：求期望（expectation）
M步：求极大（maximization)
深度在于它的数学推理涉及到比较繁杂的概率公式等，所以本文会介绍很多概率方面的知识，不懂的同学可以先去了解一些知识，当然本文也会尽可能的讲解清楚这些知识，讲的不好的地方麻烦大家评论指出，后续不断改进完善。

EM算法引入
概率模型有时候既含有观测变量，又含有隐变量或潜在变量，如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计方法估计模型参数，但是当模型含有隐变量时，就不能简单的使用这些方法，EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法，我们讨论极大似然估计，极大后验概率估计与其类似。 参考统计学习方法书中的一个例子来引入EM算法， 假设有3枚硬币，分别记做A、B、C，这些硬币正面出现的概率分别是$\pi$、$p$、$q$，进行如下实验：

先掷硬币A，根据结果选出硬币B和硬币C，正面选硬币B，反面选硬币C
通过选择出的硬币，掷硬币的结果出现正面为1，反面为0 如此独立地重复n次实验，我们当前规定n=10，则10次的结果如下所示：$$1,1,0,1,0,0,1,0,1,1$$假设只通过观测到掷硬币的结果，不能观测掷硬币的过程，问如何估计三个硬币出现正面的概率？ 我们来构建这样一个三硬币模型：$$\begin{array}{rl}P(y|\theta )& amp;=\sum _{z}P(y,z|\theta )=\sum _{z}P(z|\theta )P(y|z,\theta )\\ & amp;=\pi p^y(1-p{)}^{1-y}+(1-\pi )q^y(1-q{)}^{1-y}\end{array}$$

若$y=1$，表示这此看到的是正面，这个正面有可能是B的正面，也可能是C的正面，则$P(1|\theta )=\pi p+(1-\pi )q$

若$y=0$，则$P(0|\theta )=\pi (1-p)+(1-\pi )(1-q)$
y是观测变量，表示一次观测结果是1或0，z是隐藏变量，表示掷硬币A的结果，这个是观测不到结果的，$\theta =(\pi ,p,q)$表示模型参数，将观测数据表示为$Y=(Y_1,Y_2,...,Y_n{)}^T$，未观测的数据表示为$Z=(Z_1,Z_2,...,Z_n{)}^T$，则观测函数的似然函数是：$$\begin{array}{rl}P(Y|\theta )& amp;=\sum _{Z}P(Z|\theta )P(Y|Z,\theta )\\ & amp;=\prod _{i=0}(\pi p^{y_i}(1-p{)}^{1-y_i}+(1-\pi )q^{y_i}(1-q{)}^{1-y_i})\end{array}$$考虑求模型参数$\theta =(\pi ,p,q)$的极大似然估计，即：$$\hat{\theta }=arg\underset{\theta }{\max }logP(Y|\theta )$$这个问题没有解析解，只有通过迭代方法来求解，EM算法就是可以用于求解这个问题的一种迭代算法，下面给出EM算法的迭代过程：

首先选取初始值，记做${\theta }^0=({\pi }^0,p^0,q^0)$，第i次的迭代参数的估计值为${\theta }^i=({\pi }^i,p^i,q^i)$
E步：计算在模型参数${\pi }^i，p^i，q^i$下观测变量$y_i$来源于硬币B的概率：$${\mu }^{i+1}=\frac{{\pi }^i(p^i{)}^{y_i}(1-p^i{)}^{1-y_i}}{{\pi }^i(p^i{)}^{y_i}(1-p^i{)}^{1-y_i}+(1-{\pi }^i)(q^i{)}^{y_i}(1-p^i{)}^{1-y_i}}$$备注一下：这个公式的分母是$P(Y|\theta )$，分子表示是来源与B硬币的概率。

M步：计算模型参数的新估计值：$${\pi }^{i+1}=\frac{1}{n}\sum _{j=1}^n{\mu }_j^{i+1}$$因为B硬币A硬币出现正面的结果，所以A硬币概率就是${\mu }_j$的平均值。$$p^{i+1}=\frac{{\sum }_{j=1}^n{\mu }_j^{i+1}{y}_j}{{\sum }_{j=1}^n{\mu }_j^{i+1}}$$分子乘以$y_i$，所以其实是计算B硬币出现正面的概率。$$q^{i+1}=\frac{{\sum }_{j=1}^n(1-{\mu }_j^{i+1})y_j}{{\sum }_{j=1}^n(1-{\mu }_j^{i+1})}$$$(1-{\mu }_j^{i+1})$表示出现C硬币的概率。

闭环形成，从$P(Y|\theta )$ 到 $\pi 、p、q$一个闭环流程，接下来可以通过迭代法来做完成。针对上述例子，我们假设初始值为${\pi }^0=0.5，p^0=0.5，q^0=0.5$，因为对$y_i=1$和$y_i=0$均有${\mu }_j^1=0.5$，利用迭代公式计算得到${\pi }^1=0.5，p^1=0.6，q^1=0.6$，继续迭代得到最终的参数：$$\widehat{{\pi }^0}=0.5，\widehat{p^0}=0.6，\widehat{q^0}=0.6$$如果一开始初始值选择为：${\pi }^0=0.4，p^0=0.6，q^0=0.7$，那么得到的模型参数的极大似然估计是$\hat{\pi }=0.4064，\hat{p}=0.5368，\hat{q}=0.6432$，这说明EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。

这个例子中你只观察到了硬币抛完的结果，并不了解A硬币抛完之后，是选择了B硬币抛还是C硬币抛，这时候概率模型就存在着隐含变量！

EM算法
输入：观测变量数据Y，隐变量数据Z，联合分布$P(Y,Z|\theta )$，条件分布$P(Z|Y,\theta )$； 输出：模型参数$\theta$

(1)选择参数的初值${\theta }^0$，开始迭代
(2) E步：记${\theta }^i$为第i次迭代参数$\theta$的估计值，在第i+1次迭代的E步，计算$$\begin{array}{rl}Q(\theta ,{\theta }^i)& amp;=E_Z[logP(Y,Z|\theta )|Y,{\theta }^i]\\ & amp;=\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^i)\end{array}$$这里，$P(Z|Y,{\theta }^i)$是在给定观测数据Y和当前的参数估计${\theta }^i$下隐变量数据Z的条件概率分布；

(3) M步：求使$Q(\theta ,{\theta }^i)$极大化的$\theta$，确定第i+1次迭代的参数的估计值${\theta }^{i+1}$，$${\theta }^{i+1}=arg\underset{\theta }{\max }Q(\theta ,{\theta }^i)$$$Q(\theta ,{\theta }^i)$是EM算法的核心，称为Q函数(Q function)，这个是需要自己构造的。

(4) 重复第(2)步和第(3)步，直到收敛，收敛条件：KaTeX parse error: Expected 'EOF', got '&' at position 32: …-\theta^{i} || &̲lt; \varepsilon…或者：KaTeX parse error: Expected 'EOF', got '&' at position 58: …,\theta^{i})|| &̲lt;\varepsilon_…收敛迭代就结束了。我们来拆解一下这个M步骤，
推导逼近
主要讲解Jensen不等式，这个公式在推导和收敛都用到，主要是如下的结论：

$f(x)$是凸函数$$f(E(X))\le E(f(x))$$
$f(x)$是凹函数$$f(E(X))\ge E(f(x))$$
推导出Em算法可以近似实现对观测数据的极大似然估计的办法是找到E步骤的下界，让下届最大，通过逼近的方式实现对观测数据的最大似然估计。统计学习基础中采用的是相减方式，我们来看下具体的步骤。

增加隐藏变量$$L(\theta )=\sum _{Z}logP(Y|Z,\theta )P(Z,\theta )$$则$L(\theta )-L({\theta }^i)$为：$$\begin{array}{r}L(\theta )-L({\theta }^i)=log(\sum _{Z}P(Y|Z,{\theta }^i)\frac{P(Y|Z,\theta )P(Z,\theta )}{P(Y|Z,{\theta }^i)})-L({\theta }^i)\\ \ge \sum _{Z}P(Y|Z,{\theta }^i)log(\frac{P(Y|Z,\theta )P(Z,\theta )}{P(Y|Z,{\theta }^i)})-L({\theta }^i)\end{array}$$$\ge$这一个步骤就是采用了凹函数的Jensen不等式做转换。因为$Z$是隐藏变量，所以有KaTeX parse error: Expected 'EOF', got '&' at position 44: …P(Y|Z,\theta^i)&̲gt;0，于是继续变：
$$\begin{array}{rl}L(\theta )-L({\theta }^i)& amp;=log(\sum _{Z}P(Y|Z,{\theta }^i)\frac{P(Y|Z,\theta )P(Z,\theta )}{P(Y|Z,{\theta }^i)})-L({\theta }^i)\\ & amp;\ge \sum _{Z}P(Z|Y,{\theta }^i)log(\frac{P(Y|Z,\theta )P(Z,\theta )}{P(Z|Y,{\theta }^i)})-L({\theta }^i)\\ & amp;=\sum _{Z}P(Z|Y,{\theta }^i)log(\frac{P(Y|Z,\theta )P(Z,\theta )}{P(Z|Y,{\theta }^i)})-\sum _{Z}P(Z|Y,{\theta }^i)L({\theta }^i)\\ & amp;=\sum _{Z}P(Z|Y,{\theta }^i)log(\frac{P(Y|Z,\theta )P(Z,\theta )}{P(Z|Y,{\theta }^i)(P(Y|{\theta }^i)})\\ & amp;\ge 0\end{array}$$
也就是：$L(\theta )\ge L({\theta }^i)+{\sum }_{Z}P(Z|Y,{\theta }^i)log(\frac{P(Y|Z,\theta )P(Z,\theta )}{P(Y|Z,{\theta }^i)L({\theta }^i)})$，有下界，最大化下界，来得到近似值。这里有一个细节：$P(Y|Z,{\theta }^i)$ 变为$P(Z|Y,{\theta }^i)$？如果要满足Jensen不等式的等号，则有：$$\frac{P(Y|Z,\theta )P(Z,\theta )}{P(Y|Z,{\theta }^i)}=c$$c为一个常数，而${\sum }_{Z}P(Y|Z,{\theta }^i)=1$则：$$\begin{array}{rl}\sum _{Z}P(Y|Z,\theta )P(Z,\theta )=c\sum _{Z}P(Y|Z,{\theta }^i)& amp;=c\\ & amp;=\frac{P(Y|Z,\theta )P(Z,\theta )}{P(Y|Z,{\theta }^i)}\\ P(Y|Z,\theta )=\frac{P(Y|Z,\theta )P(Z,\theta )}{{\sum }_{Z}P(Y|Z,\theta )P(Z,\theta )}=\frac{P(Y,Z,\theta )}{P(Y,\theta )}=P(Z|Y,\theta )& \end{array}$$

大家是不是很奇怪$P(Y|Z,\theta )P(Z,\theta )$加上$\sum$之后等于什么，其实有的博客这里使用$P(Z,\theta )=P(Y^i,Z^i,{\theta }^i)$来替代$P(Y|Z,\theta )$参与计算，这样${\sum }_{Z}P(Y^i,Z^i,{\theta }^i)$，这样就方便理解来了。

于是最大化如下：

$$\begin{array}{rl}{\theta }^{i+1}& amp;=arg\underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^i)log(\frac{P(Y|Z,\theta )P(Z,\theta )}{P(Z|Y,{\theta }^i)})\\ & amp;=arg\underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^i)log(P(Y|Z,\theta )P(Z,\theta ))\\ & amp;=arg\underset{\theta }{\max }\sum _{Z}P(Z|Y,{\theta }^i)log(P(Y,Z|\theta ))\\ & amp;=arg\underset{\theta }{\max }Q(\theta ,{\theta }^i)\end{array}$$
其中$log$分母提出来是关于$Z$的${\sum }_{Z}P(Z|Y,{\theta }^i)logP(Z|Y,{\theta }^i)$，可以去掉。当然也有博客写的形式是：$$arg\underset{\theta }{\max }\sum _{i=1}^M\sum _{Z^i}P(Z^i|Y^i,{\theta }^i)log(P(Y^i,Z^i;\theta ))$$形式其实一样，表示的不一样而已。