EM算法

一、定义

EM算法，全称Expectation Maximization Algorithm，译作最大期望化算法或期望最大算法，它是一种迭代算法，从不完全数据的数据集中求解模型参数的最大似然估计方法。

其所谓“不完全数据”一般指两种情况：一种是由于观测过程本身的限制或者错误，造成观测数据成为有错漏的“不完全”数据；另一种是参数的似然函数直接优化十分困难，而引入额外的参数 (隐含的或丢失的) 后就比较容易优化，于是定义原始观测数据加上额外参数组成“完全数据”，原始观测数据自然就成为“不完全数据”。实际上，在模式识别及其相关领域, 后一种情况更为常见。

二、推导过程

（1）Jensen不等式

我们首先了解Jensen不等式（琴生不等式），它给出积分的凸函数和凹函数的积分值间的关系。那么凸凹函数的定义是什么呢？假设f是定义为实数的函数，如果对于所有的x,f(x)的二阶导数大于等于0，那么f是凸函数。当x是向量时，如果Hessian矩阵是半正定(即H>=0）,那么f是凸函数。如果f(x)的二阶导数小于0或者H<0，那么f就是凹函数。

Jensen不等于描述：

如果f是凸函数，X是随机变量，则 $E[f(X)]\geq f(E[X])$ ，当 $f''(x)\gt 0$ 或者 $H \gt 0$ 时，则 $E[f(X)]\gt f(E[X])$
如果f是凹函数，X是随机变量，则 $E[f(X)]\leq f(E[X])$ ,当 $f''(x)\lt 0$ 或者 $H \lt 0$ 时，则 $E[f(X)]\lt f(E[X])$

来看下面的图片：

Andrew Ng Lecture notes.png

f是一个凸函数，X是一个随机变量。在横轴上，X有0.5的概率取a,0.5的概率取b。X的期望是a,b值的中点（即 $\frac{a+b}{2}$ ），a,b值的期望的函数值f(E(X))。在纵轴上，f(a)和f(b)的期望是E[f(X)]。可以很直观地看出E[f(X)] > f(E[x])

（2）推导EM算法：

对于m个样本观测数据 $x =(x_1,x_2,\dots,x_m)$ 中，找出样本的模型参数 $\Theta$ ，极大化模型分布的对数似然函数
$θ=\mathop{\arg\max}_{θ} ∑_{i=1}^m log P(x_i;θ)$

如果我们得到未观察到的隐含变量 $z=(z_1,z_1,\dots,z_m)$ ，此时，我们的极大化模型参数的对数似然函数如下：
$θ=\mathop{\arg\max}_{θ} ∑_{i=1}^m log P(x_i;θ)= ∑_{i=1}^m log P(x_i,z_i;θ)(1)$

上面的例子是没有办法直接求θ 的，需要一些特殊技巧，我们首先对公式进行缩放如下：

$∑_{i=1}^m log∑_{z_i}P(x_i,z_i;θ)=∑_{i=1}^m log∑_{z(i)}Q_i(z_i)\frac{P(x(i),z(i);θ)}{Q_i(z_i)}(2)$
$≥∑_{i=1}^m∑_{z_i}Q_i(z_i)log\frac{P(x_i,z_i;θ)}{Q_i(z_i)}(3)$

上面的三个式子中，式（1）是根据联合概率密度下某个变量的边缘函数求解的，对每个样本x的所有可能z求等式右边的联合概率密度函数和，得到等式左边为随机变量X的边缘概率密度。对于（1）式直接求导非常困难，所以将每个分子分母都乘以相等的 $Q_i(z_i)$ ,转化成（2）,在等式(2)转化成等式(3)的过程，采用了Jensen不等式。

把log函数看成整体，由于Log函数的二阶导数 $-\frac{1}{x^2}< 0$ ,为凹函数。使用Jensen不等式准则(2)： $f(E[X]) \geq E[f(X)]$

我们可以把 $Q_i(z_i)$ 看成是p(x), $\frac{P(x(i),z(i);θ)}{Q_i(z_i)}$ 看成是g(x)。根据期望公式 $E[X] = \sum x*p(x)$ ,可以得到 $\sum_z Q_i(z_i)*[\frac{P(x(i),z(i);θ)}{Q_i(z_i)}]$ ，即为 $\frac{P(x(i),z(i);θ)}{Q_i(z_i)}$ 的期望

再根据凹函数对应的Jensen不等式性质
$f(E[\frac{P(x(i),z(i);θ)}{Q_i(z_i)}]) \geq E[f(\frac{P(x(i),z(i);θ)}{Q_i(z_i)}]$ ,得到公式（3）

如果要满足Jensen不等式，c为常数，则 $\frac{P(x(i),z(i);θ)}{Q_i(z_i)}=c$
由于 $Q_i(z_i)$ 是一个分布
$\sum_z Q_i(z_i) = 1$
从(2),(3)我们可以得出：
$Q_i(z_i)=\frac{P(x_i,z_i,θ)}{\sum_z P(x_i,z_i;θ)}=\frac{P(x_i,z_i,θ)}{P(x_i;θ)}=p(z_i|x_i;θ)$

等式(3)中，我们包含了隐藏数据的一个下界，如果我们能极大化这个下界，则也在尝试极大化我们对数似然函数：
$\mathop{\arg\max}_{θ}=∑_{i=1}^m∑_{z_i}Q_i(z_i)log\frac{P(x_i,z_i;θ)}{Q_i(z_i)}$
去掉上式中的常数部分，则需要极大化的对数似然下界
$\mathop{\arg\max}_{θ}=∑_{i=1}^m∑_{z_i}Q_i(z_i)log P(x_i,z_i;θ)(4)$
上式也是我们EM算法的M步， $∑_{i=1}^m∑_{z_i}Q_i(z_i)log P(x_i,z_i;θ)$ 可以看作基于条件概率分布 $Q_i(z_i)$ 的期望，是E步。

式（2）和式（3）可以写成：似然函数 $L(θ) \geq J(Z,θ)$ ,我们可以通过不断最大化J的下界，来使得 $L(θ)$ 不断提高，最终达到它的最大值。

EM极大.jpg

上图的内在含义：

首先我们固定 $θ$ ，调整 $Q(Z)$ 使下界 $J(Z,θ)$ 上升至与 $L(θ)$ 在此点处的 $θ$ 相等（黄色曲线到蓝色曲线），然后固定 $Q(Z)$ ，调整 $θ$ 使下界 $J(Z,θ)$ 达到最大值（ $θ$ 到 $θ_{t+1}$ ），然后再固定 $θ$ ,调整 $Q(Z)$ ....直到收敛到似然函数的最大值处的 $θ$

概率论中随机变量的期望计算方法

设Y是随机变量x的函数，Y=g(X)(g是连续函数）,那么
（1）X是离散性随机变量，它的分布律是 $p(X=x_k)=p_k,k=1,2,\dots$ ，若 $\sum_{k=1}^{\infty}g(x_k)p_k$ 绝对收敛，则 $E[Y]=E[g(X)]=\sum_{k=1}^{\infty} g(x_k)p_k$
（2）X是连续型随机变量，它的概率密度为f(X）,若 $\int_{-\infty}^{+\infty}g(x)f(x)dx$ 绝对收敛，则有 $E(Y) = E[g(X)]=\int_{-\infty}^{+\infty}g(x)f(x)dx$

(3)EM算法流程

输入：观测数据 $x=(x_1,x_2,\dots,x_m)$ ，联合分布 $p(x,z;θ)$ ,条件概率分布p(z|x;θ)，最大迭代次数J

随机初始化模型参数θ的初始值 $θ$
for j from 1 to j 开始EM算法迭代:
(a)E步:计算联合分布的条件概率期望：
$Q_i(z_i)= P(z_i|x_i;θ_j)$
$Q(θ,θ_j) =∑_{i=1}^m∑_{z(i)}Q_i(z_i) log P(x_i,z_i;θ_j)$
(b)M步：极大化 $L(θ,θ_j)$ ,得到 $θ_j$
$θ_{j+1}=\mathop{\arg\max}_{θ}L(θ,θ_j)$
(c)如果 $θ_{j+1}$ 已收敛，则算法结束。否则继续回步骤(a)进行E步迭代。
输出：模型参数θ

(4)EM算法收敛

假设 $P(X|\Theta)$ 为观测数据的似然函数， $\Theta^{(i)}(i=1,2,\dots)$ 为EM算法得到的参数估计序列， $p(X|\Theta^{(i)})(i=1,2,\dots)$ 为对应的参数序列，则 $p(X|\Theta^{(i)})$ 是单调递增的,即
$\sum_{i=1}^mp(x_i|θ_{i+1}) \geq \sum_{i=1}^mp(x_i|θ_{i})$

证明
$p(x_i|θ_i)=\frac{p(x_i,z_i|θ_i)}{p(z_i|x_i,θ_i)}$
取对数
$log\space p(x_i|θ_i) = log \space p(x_i,z_i|θ_i) - log\space p(z_i|x_i,θ_i)$
由于式(4)
$Q(θ,θ_i)=∑_{i=1}^m∑_{z_i}Q_i(z_i)log P(x_i,z_i;θ)=∑_{i=1}^m∑_{z_i}logP(z_i|x_i;θ_i) P(x_i,z_i;θ)$
令 $H(θ,θ_i) = ∑_{i=1}^m∑_{z_i}P(z_i|x_i;θ)log P(z_i|z_i;θ_i)$
则：
$log\space p(x_i|θ_i) = [Q(θ,θ_{i+1})-Q(θ,θ_i)] - [H(θ,θ_{i+1}),H(θ,θ_i)]$
由于 $θ_{i+1}$ 使得 $Q(θ,θ_i)$ 极大，因此有 $Q(θ,θ_{i+1})-Q(θ,θ_i) \geq 0$
而对于第二部分，我们有:
$H(θ,θ_{i+1}),H(θ,θ_i) =∑_{i=1}^m∑_{z_i}P(z_i|x_i;θ_i)log \frac{ P(z_i|z_i;θ_{i+1})}{P(z_i|z_i;θ_i)}\space(5)$
$\leq ∑_{i=1}^m log(∑_{z_i}P(z_i|x_i;θ_i) \frac{ P(z_i|z_i;θ_{i+1})}{P(z_i|z_i;θ_i)} \space(6)$
$=\sum_{i=1}^m log (\sum_{z_i} p(z_i|x_i;θ_{i+1}))=0 \space(7)$