进击数据挖掘十大算法(六):EM算法
引言
概率模型有时既含有观测变量，又含有隐变量或潜在变量。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计法，或贝叶斯估计法估计模型参数。但是 ，当模型含有隐变量时，就不能简单地使用这些估计方法。EM算法就是含有隐变量的概率模型参数的极大似然估计法，或极大后验概率估计法。
EM算法作为一种数据添加算法，在近几十年得到迅速的发展，主要源于当前科学研究以及各方面实际应用中数据量越来越大的情况下，经常存在数据缺失或者不可用的的问题，这时候直接处理数据比较困难，而数据添加办法有很多种，常用的有神经网络拟合、添补法、卡尔曼滤波法等等，但是EM算法之所以能迅速普及主要源于它算法简单，稳定上升的步骤能非常可靠地找到“最优的收敛值”。随着理论的发展，EM算法己经不单单用在处理缺失数据的问题，运用这种思想，它所能处理的问题更加广泛。有时候缺失数据并非是真的缺少了，而是为了简化问题而采取的策略，这时EM算法被称为数据添加技术，所添加的数据通常被称为“潜在数据”，复杂的问题通过引入恰当的潜在数据，能够有效地解决我们的问题。
介绍EM算法之前，我们需要介绍极大似然估计以及Jensen不等式。
一、极大似然估计
(1) 学生身高问题 我们需要调查学校的男生和女生的身高分布。 假设在校园里随机找了100个男生和100个女生。他们共200个人。将他们按照性别划分为两组，然后先统计抽样得到的100个男生的身高。假设他们的身高是服从正态分布的，但是这个分布的均值 \(\mu\) 和方差 \(\sigma^2\) 我们不知道，这两个参数就是我们要估计的。记作 \(\theta = [\mu, \sigma]^T\)。
问题：我们知道样本所服从的概率分布模型和一些样本，需要求解该模型的参数。
我们已知的有两个：样本服从的分布模型、随机抽取的样本；我们未知的有一个：模型的参数。根据已知条件，通过极大似然估计，求出未知参数。总的来说：极大似然估计就是用来估计模型参数的统计学方法。
(2) 如何估计
设样本集 \(X=x_1,x_2,\cdots,x_N\)，其中 \(N=100\) , \(p(x_i|\theta)\) 为概率密度函数，表示抽到男生的身高 \(x_i\)的概率。由于100个样本之间独立同分布，所以同时抽到这100个男生的概率就是他们各自概率的乘机，也就是样本集 \(X\) 中各样本的联合概率，用下式表示：
\[L(\theta) = L(x_1,\cdots,x_n;\theta)=\prod_{i=1}^np(x_i;\theta)\]
这个概率反映了，在概率密度函数的参数是 \(\theta\) 时，得到 \(X\) 这组样本的概率。 我们需要找到一个参数 \(\theta\) ，使得抽到 \(X\) 这组样本的概率最大，也就是说需要其对应的似然函数 \(L(\theta)\) 最大。满足条件的 \(\theta\) 叫做 \(\theta\) 的最大似然估计量，记为 ：
\[\hat{\theta} = argmax L(\theta)\]
(3) 求最大似然函数估计值的一般步骤
(a) 首先，写出似然函数：
\[L(\theta) = L(x_1,\cdots,x_n;\theta)=\prod_{i=1}^np(x_i;\theta)\]
(b) 然后，对似然函数取对数：
\[H(\theta) = \ln L(\theta) = \ln\prod_{i=1}^np(x_i;\theta) = \sum_{i=1}^n\ln p(x_i;\theta) \]
- (c) 接着，对上式待求参数求导，令导数为0，得到似然方程
- (d) 最后，求解似然方程，得到的参数 \(\theta\) 即为所求。
二、Jensen不等式
设 \(f(x)\) 是定义域为实数的函数，如果对于所有的实数 \(x\)， \(f(x)\) 的二次导数大于等于0，那么 \(f(x)\) 是凸函数(下凹)。
\(Jensen\) 不等式表述如下：如果 \(f(x)\) 是凸函数，\(X\) 是随机变量，那么：\(E[f(x)] \geq f[E(x)]\) 。当且仅当\(X\) 是常量时，上式取等号。其中，\(E(x)\) 表示 \(X\) 的数学期望。
图示中，实线 \(f(x)\) 是凸函数，\(X\) 是随机变量，有0.5的概率是 \(a\)，有0.5的概率是 \(b\)。\(X\) 的期望值就是 \(a\) 和 \(b\) 的中值了，图中可以看到 \(E[f(x)] \geq f[E(x)]\) 成立。
注：
- Jensen不等式应用于凹函数时，不等号方向反向。当且仅当 \(X\) 是常量时，Jensen不等式等号成立。
- 关于凸函数，百度百科中是这样解释的——“ 对于实数集上的凸函数，一般的判别方法是求它的二阶导数，如果其二阶导数在区间上非负，就称为凸函数（向下凸）”。关于函数的凹凸性，中国数学界关于函数凹凸性定义和国外很多定义是反的。国内教材中的凹凸，是指曲线，而不是指函数，图像的凹凸与直观感受一致，却与函数的凹凸性相反。只要记住“函数的凹凸性与曲线的凹凸性相反”就不会把概念搞乱了”。关于凹凸性这里，确实解释不统一。这里暂时以函数的二阶导数大于零定义凸函数，此处不会过多影响EM算法的理解，只要能够确定何时 \(E[f(x)] \geq f[E(x)]\) 或者 \(E[f(x)] \leq f[E(x)]\) 就可以。
三、EM算法
3.1 问题引进
我们目前有100个男生和100个女生的身高，共200个数据，但是我们不知道这200个数据中哪个是男生的身高，哪个是女生的身高。假设男生、女生的身高分别服从正态分布，则每个样本是从哪个分布抽取的，我们目前是不知道的。这个时候，对于每一个样本，就有两个方面需要猜测或者估计：
这个身高数据是来自于男生还是来自于女生？男生、女生身高的正态分布的参数分别是多少？EM算法要解决的问题正是这两个问题。
3.2 EM算法的导出
我们面对一个含有隐变量的概率模型（如上面的男女，即各自的正态分布），目标是极大化观测数据 \(Y\) 关于参数 \(\theta\) 的对数似然函数，即极大化
\[L(\theta) = \log P(Y|\theta)=\log\mathop{\sum}\limits_{Z}P(Y,Z|\theta)=\log\left(\mathop{\sum}\limits_{Z }P(Y|Z,\theta)P(Z|\theta)\right)\]
注意到这一极大化的主要困难是式中有未观测数据并有包含和的对数。
事实上，EM算法是通过迭代逐步近似极大化 \(L(\theta)\) 的。假设在第 \(i\) 次迭代后 \(\theta\) 的估计值是 \(\theta^{(i)}\)。我们希望新估计值 \(\theta\) 能使 $L() $ 增加，即 \(L(\theta)>L(\theta^{(i)})\) ，并逐步达到极大值。为此，考虑两者的差：
\[L(\theta)-L(\theta^{(i)})= \log\left(\mathop{\sum}\limits_{Z }P(Y|Z,\theta)P(Z|\theta)\right)-\log P(Y|\theta^{(i)})\]
利用Jensen不等式得到其下界：
\[L(\theta)-L(\theta^{(i)})= \log\left(\mathop{\sum}\limits_{Z}P(Z|Y,\theta^{(i)})\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})}\right)-\log P(Y|\theta^{(i)})\]
\[\geq \mathop{\sum}\limits_{Z}P(Z|Y,\theta^{(i)})\log \frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})}-\log P(Y|\theta^{(i)})\]
\[= \mathop{\sum}\limits_{Z}P(Z|Y,\theta^{(i)})\log \frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}\]
令
\[B(\theta,\theta^{(i)})\hat{=}L(\theta^{(i)})+\mathop{\sum}\limits_{Z}P(Z|Y,\theta^{(i)})\log \frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}\]
则
\[L(\theta) \geq B(\theta,\theta^{(i)})\]
即函数 \(B(\theta,\theta^{(i)})\) 是 \(L(\theta)\) 的一个下界，而且由 \(B(\theta,\theta^{(i)})\) 定义式可知，
\[L(\theta^{(i)})=B(\theta^{(i)},\theta^{(i)})\]
因此，任何可以使 \(B(\theta^{(i)},\theta^{(i)})\) 增大的 \(\theta\) ，也可以使 \(L(\theta)\) 增大。为了使 \(L(\theta)\) 有尽可能大的增长，选择 \(\theta^{(i+1)}\) 使 \(B(\theta,\theta^{(i)})\) 达到极大，即
\[\theta^{(i+1)}=arg \mathop{max}\limits_{\theta}B(\theta,\theta^{(i)})\]
现在求 \(\theta^{(i+1)}\) 的表达式，省去对 \(\theta\) 的极大化而言是常数的项，有
\[\theta^{(i+1)}=arg \mathop{max}\limits_{\theta}\left(L(\theta^{(i)})+\mathop{\sum}\limits_{Z}P(Z|Y,\theta^{(i)})\log \frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}\right)\]
\[=arg \mathop{max}\limits_{\theta}\left(\mathop{\sum}\limits_{Z}P(Z|Y,\theta^{(i)})\log P(Y|Z,\theta)P(Z|\theta)\right)\]
\[=arg \mathop{max}\limits_{\theta}\left(\mathop{\sum}\limits_{Z}P(Z|Y,\theta^{(i)})\log P(Y,Z|\theta)\right)\]
\[=arg \mathop{max}\limits_{\theta}Q(\theta,\theta^{(i)})\]
上式等价于EM算法的一次迭代，即求Q函数及其极大化。EM算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。
3.3 EM算法
- 输入：观测变量数据 \(Y\)，隐变量数据 \(Z\)，联合分布 \(P(Y,Z|\theta)\) ，条件分布 \(P(Z|Y,\theta)\)；
- 输出：模型参数 \(\theta\)
(1) 选择参数的初值 \(\theta^{(0)}\)，开始迭代；
(2) E步：记 \(\theta^{(i)}\) 为第 \(i\) 次迭代参数 \(\theta\) 的估计值，在第 \(i+1\) 次迭代的E步，计算
\[Q(\theta,\theta^{(i)})=E_z[\log P(Y,Z|\theta)|Y,\theta^{(i)}]\]
\[=\mathop{\sum}\limits_{Z}\log P(Y,Z|\theta)P(Z|Y,\theta^{(i)})\]
这里，\(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 \mathop{max}\limits_{\theta}Q(\theta,\theta^{(i)})\]
(4) 重复第(2)步和第(3)步，直至收敛
上式的函数 \(Q(\theta,\theta^{(i)})\) 是EM算法的核心，称为 \(Q\) 函数。
下面关于EM算法作几点说明：
步骤(1)：参数的初值可以任意选择，但需注意EM算法对初值是敏感的。
步骤(2)：每次迭代实际在求Q函数及其极大。
步骤(3)：M步每次迭代使似然函数增大或达到局部极值。
步骤(4)：给出停止迭代的条件，一般是对较小的正数 \(\epsilon_1,\epsilon_2\)，若满足
\(||\theta^{(i+1)}-\theta^{(i)}||<\epsilon_1\) 或 \(||Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})||<\epsilon_2\)
则停止迭代。
3.4 EM算法的收敛性
问题：EM算法是否收敛？如果收敛，是否收敛到全局最大值或局部极大值？
定理1：设 \(P(Y|\theta)\) 为观测数据的似然函数， \(\theta^{(i)}\) 为EM算法得到的参数估计序列， \(P(Y|\theta^{(i)})\) 为对应的似然函数序列，则 $ P(Y|^{(i)})$ 是单调递增的，即
\[P(Y|\theta^{(i+1)})\geq P(Y|\theta^{(i)})\]
证明：由于
\[P(Y|\theta) = \frac{P(Y,Z|\theta)}{P(Z|Y,\theta)}\]
取对数有
\[\log P(Y|\theta) = \log P(Y,Z|\theta)-\log P(Z|Y,\theta)\]
又因为
\[Q(\theta,\theta^{(i)})=\mathop{\sum}\limits_{Z}\log P(Y,Z|\theta)P(Z|Y,\theta^{(i)})\]
令
\[H(\theta,\theta^{(i)})=\mathop{\sum}\limits_{Z}\log P(Z|Y,\theta)P(Z|Y,\theta^{(i)})\]
于是对数似然函数可以写成
\[\log P(Y|\theta)=Q(\theta,\theta^{(i)})-H(\theta,\theta^{(i)})\]
在上式分别取 \(\theta\) 为 \(\theta^{(i)}\) 和 \(\theta^{(i+1)}\) 并相减，有
\[\log P(Y|\theta^{(i+1)}) - \log P(Y|\theta^{(i)})\]
\[=[Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})]-[H(\theta^{(i+1)},\theta^{(i)})-H(\theta^{(i)},\theta^{(i)})]\]
为证明定理成立，只需证明上式是非负的。上式右端第1项，由于 \(\theta^{(i+1)}\) 使 \(Q(\theta,\theta^{(i)})\) 达到极大，所以有
\[Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})\geq 0\]
其第2项，由H定义式可得：
\[H(\theta^{(i+1)},\theta^{(i)})-H(\theta^{(i)},\theta^{(i)})=\mathop{\sum}\limits_{Z}\left(\log \frac{P(Z|Y,\theta^{(i+1)})}{P(Z|Y,\theta^{(i)})}\right)P(Z|Y,\theta^{(i)})\]
\[\leq \log\left(\mathop{\sum}\limits_{Z}\frac{P(Z|Y,\theta^{(i+1)})}{P(Z|Y,\theta^{(i)})}P(Z|Y,\theta^{(i)})\right)\]
\[=\log \left(\mathop{\sum}\limits_{Z}P(Z|Y,\theta^{(i+1)})\right)=0\]
这里的不等号由Jensen不等式得到。
综上所述，右端是非负的，定理成立。
定理2：设 \(L(\theta)=\log P(Y|\theta)\) 为观测数据的对数似然函数， \(\theta^{(i)}\) 为EM算法得到的参数估计序列， \(L(\theta^{(i)})\) 为对应的对数似然函数序列。
(1) 如果 $ P(Y|)$ 有上界，则\(L(\theta^{(i)})=\log P(Y|\theta^{(i)})\) 收敛到某一值 \(L^*\) ；
(2) 在函数 \(Q(\theta,\theta^{'})\) 与 \(L(\theta)\) 满足一定条件下，由EM算法得到的参数估计序列 \(\theta^{(i)}\) 的收敛值 \(\theta^*\) 是 \(L(\theta)\) 的稳定点。
证明：
- \(L(\theta)=\log P(Y|\theta)\) 的单调性及 $ P(Y|)$ 的有界性
- 证明从略
四、Reference
- 李航《统计学习方法》