EM算法

本文介绍了机器学习领域经典算法之一的EM算法,也被称为数据挖掘十大算法之一。EM算法是一种迭代算法,用于含有隐变量的概率模型参数的极大似然估计或极大后验概率估计。EM算法的每次迭代由两步组成:E步,求期望(expectation);M步,求极大(maximization)。所以这一算法称为期望极大算法(expectation maximization algorithm),简称为EM算法。


预知识

符号说明

符号 含义
$D$ 已有的数据(data)
$\theta$ 要估计的参数(parameter)
$p(\theta)$ 先验概率(prior)
$p(\theta\mid{D})$ 后验概率(posterior)
$p(D)$ 数据分布(evidence)
$p(D\mid\theta)$ 似然函数(likelihood of θ w.r.t. D)
$p(y,\theta\mid{D})$ 已知数据条件下的y,θ概率

先验概率&后验概率

先验概率(prior)与后验概率(posterior)简称为先验后验

这两个概念来自于贝叶斯定理,我们肯定是针对同一个事物才有先后之分,如果针对两个事物,先后不是没有意义了么?那这个共同的对象,就是我们的参数$\theta$。后验概率是指掌握了一定量的数据后我们的参数分布是怎么样的,表示为$p(\theta|D)$;那先验就是在没有掌握数据时我们的参数怎么分布。

看到这里,你可能会问:如果连数据都没有,我怎么知道我的参数是怎么分布的?你提出这个问题,就说明你是一个赤裸裸的频率派学家,你需要通过数据来得到你的参数!而这并不是贝叶斯派的考虑,贝叶斯估计最重要的就是那个先验的获得。虽然你这次的一组数据,比如说扔三次硬币产生的序列是(110),但是其实根据我历史的经验来看,一枚硬币正反面其实很有可能是按照均匀分布来的,只不过可能因为你抛得次数少了所以产生了不是均匀分布的效果,因此我要考虑以往的经验在里面。

你可能又会问:那你这个均匀分布不就是完全猜来的嘛,你怎么知道我这次是不是一样的硬币呢?没错!就是“猜来的”。先验在很多时候完全是假设,然后去验证有的数据是否吻合先验猜想,所以这里的猜很重要。还要注意,先验一定是与数据无关的,你不能看到了数据再做这些猜想,一定是没有任何数据之前你就猜了一个参数的先验概率。

贝叶斯定理

通常,事件A在事件B(发生)的条件下的概率,与事件B在事件A(发生)的条件下的概率是不一样的;然而,这两者是有确定的关系的,贝叶斯定理就是这种关系的陈述。贝叶斯公式的一个用途在于通过已知的三个概率函数推出第四个。

作为一个普遍的原理,贝叶斯定理对于所有概率的解释是有效的;然而,频率主义者贝叶斯主义者对于在应用中,某个随机事件的概率该如何被赋值,有着不同的看法:频率主义者根据随机事件发生的频率,或者总体样本里面的发生的个数来赋值概率;贝叶斯主义者则根据未知的命题来赋值概率。这样的理念导致贝叶斯主义者有更多的机会使用贝叶斯定理。
$$
P(A|B) = \frac{P(B|A)P(A)}{P(B)}\tag{1}
$$
在贝叶斯定理中,每个名词都有约定俗成的名称:

  • P(A|B)是已知B发生后A的条件概率,也因为和B有关被称作A的后验概率。
  • P(B|A)是已知A发生后B的条件概率,也因为和A有关被称作B的后验概率。
  • P(A)是A的先验概率(或边缘概率),因为它不考虑任何B方面的因素。
  • P(B)是B的先验概率(或边缘概率),因为它不考虑任何A方面的因素。

按这些术语,贝叶斯定理可表述为:
$$
后验概率 = \frac{相似度\times先验概率}{标准化向量}\tag{2}
$$

也就是说,后验概率与先验概率和相似度的乘积成正比。

极大似然估计

极大似然估计法认为参数是固有的,但是可能由于一些外界噪声的干扰,使数据看起来不是完全由参数决定的。没关系,数学家们觉得,虽然有误差存在,但只要让在这个数据给定的情况下,找到一个概率最大的参数就可以了。那问题其实就变成了一个条件概率最大的求解,即求使得$p(\theta|D)$ 最大的参数$\theta$,形式化表达为求解
$$
\arg\max_{\theta}p(\theta|D)\tag{3}
$$

而根据条件概率公式有
$$
p(\theta|D)=\frac{p(D|\theta)p(\theta)}{p(D)}\tag{4}
$$

因为我们在极大似然估计中假设是确定的,所以$p(\theta)$就是一个常数。$p(D)$同样是根据已有的数据得到的,也是确定的,或者我们可以把其看作是对整个概率的一个归一化因子。这时候,求解公式(3)就变成了求解
$$
\arg\max_{\theta}p(D|\theta)\tag{5}
$$
(5)式中的$p(D|\theta)$就是似然函数,我们要做的就是求一个是似然最大的参数,所以称为极大似然估计。想求解这个问题,需要假设我们的数据是相互独立的。$D=\{y1,y_2,y_3,…,y_n\}$,这时候有
$$
p(D|\theta)=\prod_{i=1}^{n}p(y_i|\theta)\tag{6}
$$
一般对(6)式取对数求解对数极大似然,就可以把连乘变成求和,然后求导取极值点就是要求的参数值,不在此赘述。

Jensen不等式

若f(x)是区间(a,b)上的凹函数,则对任意的$x_1,x_2,x_3,…,x_n\in(a,b)$,有:
$$
f(\frac{x_1+x_2+x_3+…+x_n}{n})\le \frac{f(x_1)+f(x_2)+f(x_3)+…+f(x_n)}{n}\tag{7}
$$
当且仅当$x_1=x_2=x3=…=x_n$时等号成立。

若f(x)是区间(a,b)上的凸函数,则对任意的$x_1,x_2,x_3,…,x_n\in(a,b)$,有:
$$
f(\frac{x_1+x_2+x_3+…+x_n}{n})\ge \frac{f(x_1)+f(x_2)+f(x_3)+…+f(x_n)}{n}\tag{8}
$$
当且仅当$x_1=x_2=x3=…=x_n$时等号成立。

其加权形式为,若f(x)是区间(a,b)上的凹函数,则对任意的$x_1,x_2,x_3,…,x_n\in(a,b)$, $\sum_{i=1}^na_i=1$有:
$$
f(a_1x_1+a_2x_2+a_3x_3+…+a_nx_n)\leq a_1f(x_1)+a_2f(x_2)+a_3f(x_3)+…+a_nf(x_n)\tag{8}
$$
凸函数同理。


EM算法

引入

概率模型有时既含有观测变量,又含有隐变量或潜在变量,如果概率模型的变量都是观测变量,那么给定数据,可以直接用极大似然估计法,或贝叶斯估计法估计模型参数。但是,当模型含有隐变量时,就不能简单的使用这些估计方法。EM算法就是含有隐变量的概率模型参数的极大似然估计法,或极大后验概率估计法。这里仅讨论极大似然估计法。

在理解EM算法之前,先来举个例子,我们叫它“三硬币模型”:现在有A、B、C三个不同的硬币,它们的材料并非均匀,投掷一次正面朝上的概率分别为$\pi、p、q$。在每轮开始前首先掷硬币A,根据结果选择另外两个硬币,如果是正面则选择B来掷,如果是反面则选择C来掷。选择完之后开始掷选择的硬币,如果是正面则记为1,反面则记为0。在每轮开始前都先掷A进行选择,然后记录B或C的投掷结果。独立重复n次试验,这里n=10,得到结果如下:
$$
1,1,0,1,0,0,1,0,1,1
$$
假设我们只能观测到掷硬币的结果,不能观测掷硬币的过程,问如何估计三硬币出现的概率,即三硬币模型的参数。

这里的观测变量Y即如上的观测结果,隐变量Z就是每次掷硬币选择的是B还是C,我们无法得知Z的值,因为我们无法观测掷硬币的过程。这种情况下,我们就可以用EM算法来解决这个问题。

感性的理解

这部分有一篇论文个人感觉写的特别好,强烈推荐大家:

What is the expectation maximization algorithm?

Do, C. B., & Batzoglou, S. (2008). What is the expectation maximization algorithm?. Nature biotechnology, 26(8), 897.

论文中提到的例子,我用python做了简单的实现,代码放在了我的GitHub上,地址:

https://github.com/zangbo/MachineLearning/tree/master/EM

关于论文的中文解析可以参考这篇博文:

如何感性的理解EM算法? - 简书

数学推导

这里我们令观测数据为Y。我们面对一个概率模型,用极大似然估计法来做的话就是要极大化观测数据Y关于参数$\theta$的对数似然函数$L(\theta)=logP(Y|\theta)$。然而在存在隐变量Z的情况下,完整数据除了Y之外还包括Z,这时我们把上面的对数似然函数就变成了:
$$
L(\theta)=logP(Y|\theta)=log\sum_ZP(Y,Z|\theta)\tag{9}
$$
极大化该似然函数是比较困难的,因为其中包含了隐变量Z,属于未观测数据。这时候我们想,我们能不能找出来一个式子,使得$L(\theta)$总能大于或等于该式子,那么如果我们能求出该式子的极大值,是不是就可以进而得出$L(\theta)$的极大值呢?我们首先对如上式子做一个等价变换:
$$
log\sum_ZP(Y,Z|\theta)=log\sum_ZG(Z)\frac{P(Y,Z|\theta)}{G(Z)}\tag{10}
$$
这里的G函数我们可以认为是只和Z有关的概率函数。接下来要用到Jensen不等式,上面已经介绍过Jensen不等式,这里我们用的是$log\sum_j\lambda_jy_j\ge\sum_j\lambda_jlogy_j$,其中$\lambda_j\ge0$,$\sum_j\lambda_j=1$。根据Jensen不等式,(10)式可以写成:
$$
log\sum_ZG(Z)\frac{P(Y,Z|\theta)}{G(Z)}\ge \sum_ZG(Z)log\frac{P(Y,Z|\theta)}{G(Z)} \tag{11}
$$
这里G需要满足的条件是$\sum_ZG(Z)=1$,接着我们想等号成立的条件是什么,也就是左式的下限是什么。使得等号成立,则必须满足$\frac{P(Y,Z|\theta)}{G(Z)}=c$,其中c为常数。那么同时满足这两个条件的G,我们可以构造它为:
$$
\begin{split}
G(Z)=&\frac{P(Y,Z|\theta^{(i)})}{\sum_ZP(Y,Z|\theta^{(i)})}\\
=&\frac{P(Y,Z|\theta)}{P(Y|\theta^{(i)})}\\
=&P(Z|Y,\theta^{(i)})
\end{split}
\tag{12}
$$
G(Z)的数学含义可以看出是在给定观测数据Y和当前参数估计$\theta^{(i)}$下隐变量数据Z的条件概率分布,代入(11)式得:
$$
\sum_ZG(Z)log\frac{P(Y,Z|\theta)}{G(Z)}=\sum_ZP(Z|Y,\theta^{(i)})\log\frac{P(Y,Z|\theta)}{P(Z|Y,\theta^{(i)})} \tag{13}
$$
此时只需要极大化右式即可,省去对$\theta$的极大化而言是常数的项,最终转化为求下式:
$$
\begin{split}
&\arg\mathop{\max} \limits_{\theta}\lgroup\sum_ZP(Z|Y,\theta^{(i)})\log\frac{P(Y,Z|\theta)}{P(Z|Y,\theta^{(i)})} \rgroup \\
=&\arg\mathop{\max} \limits_{\theta}\lgroup \sum_ZP(Z|Y,\theta^{(i)})logP(Y,Z|\theta) \rgroup \\
=&\arg\mathop{\max} \limits_{\theta}\lgroup \sum_ZG(Z)logP(Y,Z|\theta) \rgroup\\
=&\arg\mathop{\max} \limits_{\theta}\lgroup Q(\theta,\theta^{(i)}) \rgroup
\end{split}
\tag{14}
$$
可以看出,$Q(\theta,\theta^{(i)})=E_Z[logP(Y,Z|\theta)|Y,\theta^{(i)}]$,即最终转化为求完全数据的对数似然函数$logP(Y,Z|\theta)$关于在给定观测数据Y和当前参数$\theta^{(i)}$下对未观测数据Z的条件概率分布$P(Z|Y,\theta^{(i)})$的期望并且求其极大化。

这就完成了EM算法的一次迭代。由此可以得出EM算法的步骤:

  1. 选择参数的初值$\theta^{(0)}$,开始迭代;

  2. E步:记$\theta^{(i)}$为第i次迭代参数$\theta$的估计值,在第i+1次迭代的E步,计算$G(Z)=P(Z|Y,\theta^{(i)})$,并代入(14)式得到
    $$
    Q(\theta,\theta^{(i)})=\sum_ZG(Z)logP(Y,Z|\theta)\tag{15}
    $$

  3. M步:求使$Q(\theta,\theta^{(i)})$极大化的$\theta$,确定第i+1次迭代的参数估计值$\theta^{(i+1)}$
    $$
    \theta^{(i+1)}=\arg\mathop{\max} \limits_{\theta}\lgroup Q(\theta,\theta^{(i)}) \rgroup\tag{16}
    $$

  4. 重复第2步和第3步,直到收敛。

由此可以看出,EM算法是通过不断求解下界的极大化逼近求解对数似然函数极大化的算法。

下图可以看出EM算法的直观解释:

从图可以推断出EM算法不能保证找到全局最优值。事实上,它只能保证找到稳定点。因此在应用中,初值的选择变得非常重要,常用的方法是选择几个不同的初值进行迭代,然后对得到的各个估值加以比较,从中选择最好的。


参考