期望最大化算法(expectation maximization,EM)

本文主要是介绍期望最大化算法(expectation maximization,EM)，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

EM算法是一种迭代算法，用于含有隐变量(hidden variable)或缺失数据(incomplete-data)的概率模型参数的极大似然估计。EM算法被广泛用于高斯混合模型(Guassian Mixture Model,GMM)和隐马尔可夫模型(Hidden Markov Model,HMM)的参数估计。每次迭代分为两步：E步，求期望；M步，求极大。概率模型有时既包含观测变量(observable variable)有包含隐变量。

作为机器学习十大算法之一，EM算法肯定有两把刷子。EM算法是基于极大似然估计(Maximum Likelihood Estimation,MLE)理论的优化算法，既然如此，在学习EM算法之前首先要对似然函数(Likelihood function)有个了解，对于似然函数，我也一直不是特别理解，通过查阅资料，对似然函数有了比较清晰的认识，既然EM算法是极大似然估计的一种扩展，那么我们先从极大似然估计说起。

1 极大似然估计

1.1似然与概率

1.2 例子

2. EM算法

2.1 高斯混合模型(Gaussian misture model,GMM)

2.2 EM算法的原理

2.2.1 EM 算法的收敛性

3. EM算法的应用

参考文献

1 极大似然估计

引用维基百科里面的话说：似然函数是一种关于统计模型中的参数的函数，表示模型参数中的似然性。还是不太理解？只看这句话我也是云里雾里。我们来一点点的看。

1.1似然与概率

似然与概率有着千丝万缕的联系，在英语语境中，likelihood和probability是可以互换的，但是在数理统计中，这两个概念是不同的。“xxx的概率”中的xxx只能是概率空间中的事件；而“xxx的似然”中的xxx只能是参数。也就是说概率是随机变量的函数，而似然函数是参数的函数。

似然可以理解为条件概率的逆反。已知某个参数 $\large \theta$ 时，事件X会发生的概率写作： $\large P(X|\theta) = \frac{P(X,\theta)}{P(\theta)}$

由贝叶斯公式得： $\large P(\theta|X) = \frac{P(X|\theta)P(\theta)}{P(X)}$ 因此，我们可以反过来构造表示似然性的方法：已知有事件X发生，运用似然函数 $\large L(\theta|X)$ ，我们估计参数 $\large \theta$ 的可能性。在形式上，似然函数也是条件概率函数，但是我们关注的变量改变了。在贝叶斯公式中， $\large P(\theta|X)$ 表示后验概率， $\large P(X|\theta)$ 表示似然， $\large P(\theta)$ 表示先验概率，而 $\large P(X)$ 是用于归一化的证据因子，给定X后，可以把它看做常量。

1.2 例子

以抛硬币为例，抛10次硬币，实验结果中有6次为“国徽”。我们假设以 $\large \theta = 0.5$ 为参数，实验结果的概率为 $\large C_{10}^{6}0.5^{6}(1-0.5)^{4}\approx0.21$ ；我们再假设以 $\large \theta = 0.6$ 为参数，实验结果的概率为： $\large C_{10}^{6}0.6^{6}(1-0.6)^{4}\approx0.25$ ，可以看到以 $\large \theta = 0.6$ 作为参数时，抛10次硬币出现6次“国徽”的概率更大。

假设硬币的参数为 $\large \theta$ ，可以得到似然函数为： $\large L(\theta) = C_{10}^{6}\theta^{6}(1-\theta)^{4}$ ,图像如下：

可以看到似然函数在0.6处的取值最大，实际上似然函数就是推测的参数的分布。极大似然估计(Maximum likelihood estimation，MLE)是在做什么样的一件事呢，我们在知道事件发生(这里的事件就是抛10次硬币，出现6次“国徽”)的条件下，利用似然函数求出使似然函数最大的参数 $\large \theta^{\star}$ ,换句话说，极大似然估计能够找到模型的参数，使得观测发生的概率最大，反应到这个例子中就是，极大似然估计能找到模型参数 $\large \theta^{\star}$ ，当 $\large \theta^{\star} = 0.6$ 时出现事件(抛10次硬币，出现6次“国徽”)的概率最大。

但是极大似然估计并不是能估算所有概率模型的参数，比如三硬币模型、高斯混合模型等等，因为在这些模型是由观测变量、隐变量以及参数 $\theta$ 共同决定的，接下来就以高斯混合模型为例子引入EM算法。

（三硬币模型）：假设有3枚硬币，分别记作A,B,C.这些硬币正面出现的概率分别为 $\large \pi,p$ 和 $\large q$ 。进行如下试验：先掷硬币A，根据其结果选出硬币B或硬币C，正面选硬币B，反面选硬币C；然后掷选出的硬币，掷硬币的结果，出现正面记作1，出现反面记作0；独立的重复n次试验(此处n = 10)，观测结果如下： {1,1,0,1,0,0,1,0,1,1}。假设只能观测到掷硬币的结果，不能观测掷硬币的过程。那么该如何估计三硬币正面出现的概率，即三硬币模型的参数。

2. EM算法

2.1 高斯混合模型(Gaussian misture model,GMM)

高斯分布也称作正态分布，是我们生活中最常见的概率分布。我们的身高在样本数足够多的情况下就是服从高斯分布的，假设我们统计500个人(随机选取)的身高，并将身高数据绘制成直方图，如下图所示。

我们可以认为身高数据是由上面这样的高斯分布随机生成的，如果考虑到男女，那么身高数据绘制成的直方图会变成下面这种形式：

这时候，身高数据是由两个不同的高斯模型共同生成的，也就是高斯混合模型如果随机给我们一个人身高数据，让我们判断这个人的性别，比如180cm，我们可以从上面的图中看到，180cm的数据由男性的高斯分布产生的概率更大一些，所以我们可以认为这个人是男性。这其实就是高斯混合模型可以用来聚类的基本思想。

(高斯混合模型)：高斯混合模型是指具有如下形式的概率分布模型：

$P(x|\theta) = \sum_{k = 1}^{K} \alpha_{k}N(x|\theta_{k})$ (2.1)

其中， $\alpha_{k}$ 是第k个高斯模型的系数，也可以看做是权重， $\alpha_{k} \geqslant 0$ ， $\sum_{k = 1}^{K}\alpha_{k} = 1$ ， $N(x|\theta_{k})$ 是高斯分布密度， $\theta_{k} = (\mu_{k},\sigma_{k}^{2})$ 。在男女身高这个例子里， $k = 2$ 。

假设观测数据 $(x_{1},x_{2},...,x_{N})$ 单个高斯模型生成的，可以直接用极大似然估计来计算参数：

$\hat\theta = \mathop{\arg\max}_{\theta}logP(X|\theta)\\ = \mathop{\arg\max}_{\theta}L(\theta|X)\\ = \mathop{\arg\max}_{\theta}[\sum_{i= 1}^{N}log(N(x_{i}|\theta)] \\$ (2.2)

其中 $\theta = (\mu,\sigma)$ ，我们可以直接对对数似然函数求偏导让偏导数为0，求得参数。

假设观测数据 $(x_{1},x_{2},...,x_{N})$ 是由高斯混合模型生成的，此时的参数为 $\theta = (\alpha_{1},\alpha_{2},...,\alpha_{K};\theta_{1},\theta_{2},...,\theta_{K})$ ,这时候，就有

$\hat\theta = \mathop{\arg\max}_{\theta}logP(X|\theta)\\ = \mathop{\arg\max}_{\theta}L(\theta|X)\\ = \mathop{\arg\max}_{\theta}[\sum_{i= 1}^{N}log(\sum_{k = 1}^{K}\alpha_{k}N(x_{i}|\theta_{k}))]$ (2.3)

此时的对数似然函数是远比单高斯模型的复杂，而且log函数里面还包含求和，这就比较复杂了，如果还是按照对对数似然函数求偏导的方法估计参数，可以实现，但是是相当困难的，方法总是比问题多，我们这是可以通过迭代的方法求解，这就类似于非概率模型中最小化平方误差函数时使用的梯度下降法，这种算法就是EM算法。

2.2 EM算法的原理

在高斯混合模型中观测数据是由高斯模型产生的，但是是由哪一个高斯模型产生的我们是不知道的，我们引入隐变量Z，对于EM算法我们有如下定义

$\theta^{(i+1)} = \mathop{\arg\max}_{\theta}\sum_{Z}logP(X,Z|\theta)P(Z|X,\theta^{(i)})$ (2.4)

当连续两次迭代 $\theta^{(i+1)}$ 和 $\theta^{(i)}$ 的差小于任意小的正数 $\varepsilon$ 时，停止迭代，即 $||\theta^{(i+1)} -\theta^{(i)}|| < \varepsilon$ 时停止迭代。

给本来就复杂的式子再另外引入隐变量Z，这能解决问题？答案是肯定的，因为隐变量可以告诉我们观测数据 $x_{i}$ 是由哪个模型产生，这样的话我就相当于把某个数据 $x_{i}$ 择(zhai)出来，这样就回到了单个高斯模型的求参问题。但是加入隐变量有一个要求就是边缘分布不能改变，用式子表示是：

$P(x_{i}) = \sum_{z_{i}}P(x_{i}|z_{i})P(z_{i})$ (2.5)

这里的 $\large z_{i}$ 是离散分布，所以用求和符号得到边缘分布，如果是连续分布，则需要积分，而高斯混合模型中的权重 $\alpha_{k}$ 就是 $\large z_{i}$ 的先验。

2.2.1 EM 算法的收敛性

“走得太远了，不要忘记我们为什么要出发”。我们训练概率模型的目的就是求使对数似然函数 $L(\theta)$ 最大的参数 $\theta^{\star}$ ,而EM算法是通过迭代逐步近似 $L(\theta)$ 的极大值的，也就是说 $L(\theta^{(i+1)}) \geq L(\theta^{(i)})$ 那么我们怎么保证每一次迭代对应的对数似然函数值都会比上一次迭代对应的对数似然函数值要大呢。接下来我们证明一下。

对数似然函数可以写成下面的式子：

$L(\theta) = logP(X|\theta) = log\frac{P(X,Z|\theta)}{P(Z|X,\theta)}\\ = logP(X,Z|\theta) - logP(Z|X,\theta)$ (2.6)

对上式等式两边分别以 $P(Z|X,\theta^{(i)})$ 为分布对函数求期望得：

$E_{P(Z|X,\theta^{(i)})}[logP(X|\theta)] = E_{P(Z|X,\theta^{(i)})}[logP(X,Z|\theta) - logP(Z|X,\theta)]$ (2.7)

等式左边为： $E_{P(Z|X,\theta^{(i)})}[logP(X|\theta)] \\=logP(X|\theta)\sum_{Z}P(Z|X,\theta^{(i)})\\ =logP(X|\theta)$ (2.8)

将式(2.8)代入(2.7)并式(2.7)等号右边展开得：

$logP(X|\theta) = \sum_{Z}logP(X,Z|\theta)P(Z|X,\theta^{(i)}) - \sum_{Z}log(Z|X,\theta)P(Z|X,\theta^{(i)})$ (2.9)

分别令:

$Q(\theta,\theta^{(i)}) = \sum_{Z}logP(X,Z|\theta)P(Z|X,\theta^{(i)})$ (2.10)

$H(\theta,\theta^{(i)}) = \sum_{Z}logP(Z|X,\theta)P(Z|X,\theta^{(i)})$ (2.11)

其中 $Q(\theta,\theta^{(i)})$ 中的第一个变量表示要极大化的参数，第二个变量表示参数的当前估计值。而且这个 $Q(\theta,\theta^{(i)})$ 就是EM算法定义中要最大化的项，那么比较两个数值大小的方法通常有两种，一种是做差，一种是做商，这里我们采用做差的方法来验证。对上式中的 $\large \theta$ 分别取 $\large \theta^{(i)}$ 和 $\large \theta^{(i+1)}$ 并相减得到：

$logP(X|\theta^{(i+1)}) - logP(X|\theta^{(i)}) \\ =[Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})]-[H(\theta^{(i+1)},\theta^{(i)})-H(\theta^{(i)},\theta^{(i)})]$ (2.12)

由式(2.4)知道，由于 $\theta^{(i+1)}$ 是使 $Q(\theta,\theta^{(i)})$ 最大化的值，所以 $Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)}) \geq 0$ 。对于第二项，有：

所以就得到了： $logP(X|\theta^{(i+1)}) - logP(X|\theta^{(i)}) \geq 0$ ，这样就证明了每一次迭代比上一次迭代得到的对数似然函数要大，这样EM算法才能收敛。

关于EM算法有几点需要说明：

1、因为EM算法是一种迭代算法，所以需要给参数设置初值，初值的选择是随意的，但是EM算法对初值是敏感的，究竟选择怎么样的初值使模型训练结果最好，这需要不断的试

2、前面说到，EM算法中，E步是求期望，那么求得是谁的期望呢。求得是联合分布 $P(X,Z|\theta)$ 的对数似然函数以 $P(Z|X,\theta^{(i)})$ 为分布的期望，即 $E_{P(Z|X,\theta^{(i)})}(logP(X,Z|\theta)) = \sum_Z logP(X,Z|\theta))P(Z|X,\theta^{(i)})$ ，把它定义为 $Q(\theta,\theta^{(i)})$ 。

3、M步就是求期望的最大值，也就是 $Q(\theta,\theta^{(i)})$ 的最大值，得到 $\theta^{(i+1)}$ ，这样就完成一次迭代。

4、当连续两次迭代 $\theta^{(i+1)}$ 和 $\theta^{(i)}$ 的差小于任意小的正数 $\varepsilon$ 时，停止迭代，即 $||\theta^{(i+1)} -\theta^{(i)}|| < \varepsilon$ 时停止迭代。

3. EM算法的应用

以三硬币模型为例，使用EM算法求解模型参数。为了方便阅读，把三硬币模型的描述搬下来。

E步：

$P(X|\theta) = \prod_{j = 1}^{n}[\pi p^{x_{j}}(1-p)^{1-x_{j}} + (1-\pi)q^{x_{j}}(1-q)^{1-y_{j}} ]$

计算在模型参数 $\pi^{(i)},p^{(i)},q^{(i)}$ 下观测数据 $x_{j}$ 来自掷硬币B的概率

$\mu^{(i+1)}_{j} = \frac{\pi^{(i)}(p^{(i)})^{y_{j}}(1-p^{(i)})^{1-y_{j}}}{\pi^{(i)}(p^{(i)})^{y_{j}}(1-p^{(i)})^{1-y_{j}} + (1-\pi^{(i)})(q^{(i)})^{y_{j}}(1-q^{(i)})^{1-y_{j}}}$