期望极大化（EM，Expectation Maximization）是一种迭代算法，1977年由Dempster等人总结提出，用于含有隐变量（hidden variable）的概率模型参数的极大似然估计或极大后验概率估计。

概率模型有时既含有观测变量（observable variable）,又含有隐变量或潜变量（latent variable）。如果概率模型的变量都是观测变量，那么给定数据后，可以直接用极大似然估计或贝叶斯估计；但当模型含有隐变量时（如一个词所属的主题，聚类中的样本列表），难以求得参数的解析解，就需要用到EM算法。

EM算的每次迭代由两步组成：

• E步：求期望步（expectation）
• M步：求极大（maximization）

1 EM算法的引入

1.1 E.g.1 两硬币模型

（1）一般情况

假设有2枚硬币A，B，正面朝上的概率分别为 $p_1,p_2$。为了估计正面朝上的概率，每次取一枚硬币，连续抛5次，记录结果如下：
则， $p_1 = 6/(6+9)=0.4$， $p_2=5/(5+5)=0.5$

（2）含有隐变量的情况

若每次使用哪枚硬币抛掷未知，结果如下：此时相当于增加了隐变量 $Z=(z_1,z_2,z_3,z_4,z_5)$，其中， $z_i=0$代表第 $i$轮使用硬币A抛掷， $z_i=1$代表使用硬币B抛掷。

为了估计出 $p_1,p_2$首先需要估计出 $Z$，由于 $p_1,p_2$未知，可先随机初始化 $p_1,p_2$，估计出 $Z$，再利用 $Z$估计新的 $p_1,p_2$。

方法一：

令 $p_1=0.2,p_2=0.7$，则在第一轮掷硬币中，

• 硬币A，3正2反的概率： $0.2 *0.2*0.2*0.8*0.8 = 0.00512$
• 硬币B，3正2反的概率： $0.7 *0.7*0.7*0.3*0.3 = 0.03087$

根据最大似然估计判定得到 $Z=(1,0,0,1,0)$，硬币A，5正10反；硬币B，6正4反。
则 $p_1=0.33,p_2=0.6$，相对于初始值更加接近真实值了。

可按照上述思路，用估计得到的 $p_1,p_2$再来估计 $Z$（E-步），再用 $Z$来估计新的 $p_1,p_2$（M-步），反复迭代下去，直到 $p_1,p_2$的值不再改变。

方法二：

上述方法中，根据每一轮使用A和B的概率，简单判定为非A即B，可不作出判定，而保留使用A和B的概率，即隐变量也是存在分布的。

如第1轮，

• 使用硬币A的概率： $0.00512/(0.00512+0.03087)=0.14$
• 使用硬币B的概率： $1-0.14=0.86$
  依次计算出其他4轮的概率如下：

此时，完成了对于隐变量 $Z$的概率分布的估计（E-步），再按照极大似然来估计 $p_1，p_2$（M-步）。

如针对硬币A，第1轮3正2反，相当于正面概率为 $0.14*3=0.42$，反面概率为 $0.14*2=0.28$，5轮列表如下：

此时， $p_1=4.22/(4.22+7.98)=0.35$

Do C B , Batzoglou S 在《What is the expectation maximization algorithm？》中提及的类似例子如下：

1.2 E.g.2 三硬币模型

假设有3枚硬币A，B，C，这些硬币正面向上的概率分别为 $\pi,p,q$。
进行如下抛硬币实验：先抛硬币A，根据其结果选择硬币B或C：正面选B，反面选C；然后抛出选择的硬币，出现正面记作1，出现反面记作2。

独立重复 $n(n=10)$次试验，观测结果如下： $1101001011$根据观测结果，估计三枚硬币正面向上的概率，即三枚硬币模型参数 $\theta = (\pi,p,q)$。
上述模型可以表示为：

$\begin{aligned} P(x \mid \theta) = & \sum_{z} P(x,z \mid \theta) \\ = &\sum_{z} P(z \mid \theta)P(x \mid z,\theta )\\ = &\pi p^x(1-p)^{1-x}+(1-\pi) q^x(1-q)^{1-x}\end{aligned}$

其中，随机变量 $x$是观测变量，表示本次试验观测结果是1或0；随机变量 $z$是隐变量，表示未观测到的抛掷硬币A的结果。
注意：随机变量 $x$的数据可以观测，随机变量 $z$的数据不可观测。

令观测数据 $X=(x_1,x_2,...,x_n)$，未观测数据 $Z=(z_1,z_2,...,z_n)$，则观测数据的似然函数：

$\begin{aligned} P(X \mid \theta) = &\sum_{Z} P(z \mid \theta)P(X \mid z,\theta) \\ = & \prod_{i=1}^{n}(\pi p^{x_i}(1-p)^{1-x_i}+(1-\pi) q^{x_i}(1-q)^{1-x_i})\end{aligned}$求参数 $\theta = (\pi,p,q)$的极大似然估计： $\hat \theta = \arg \max_\theta \log P(X \mid \theta)$

上述问题没有解析解，只能通过迭代的方法求解：选取参数的初始值，记作 $\theta_0=(\pi_0,p_0,q_0)$，不断迭代计算参数的估计值，直至收敛。第 $j$次迭代，参数估计值为 $\theta_j=(\pi_j,p_j,q_j)$。其中，第 $j+1$次迭代如下：

• E步：计算模型在参数 $\theta_j=(\pi_j,p_j,q_j)$下观测数据 $y_i$来自硬币B的概率:

$\mu_{(i,j+1)}= \frac{\pi_j p_j^{x_i}(1-p_j)^{1-x_i}}{\pi_j p_j^{x_i}(1-p_j)^{1-x_i}+(1-\pi_j) q_j^{x_i}(1-q_j)^{1-x_i}}$

• M步：计算新的参数估计值:
  $\pi_{j+1}=\frac{1}{n}\sum_{i=1}^{n}\mu_{(i,j+1)}$

$p_{j+1} = \frac{\sum_{i=1}^{n}\mu_{(i,j+1)}x_i}{\sum_{i=1}^{n}\mu_{(i,j+1)}}$

$q_{j+1} = \frac{\sum_{i=1}^{n}(1-\mu_{(i,j+1)})x_i}{\sum_{i=1}^{n}(1-\mu_{(i,j+1)})}$

假设模型参数 $\theta_0=(\pi_0,p_0,q_0)=(0.5,0.5,0.5)$，对于 $x_i=1与x_i=0，i=1,2,...,10$均有 $\mu_{(i,1)}=0.5$。
根据迭代，可得到 $\theta_1=(0.5,0.6,0.6)$则 $\mu_{(i,2)}=0.5,i=1,2,...,10$继续迭代得到， $\theta_2=(0.5,0.6,0.6)$模型收敛，得到参数 $\theta$的极大似然估计： $\hat \theta=(0.5,0.6,0.6)$ $\pi=0.5$表示硬币A均匀， $p=0.6,q=0.6$符合对于数据的直观观察，10次抛硬币中，正面出现了6次。
如果初始值 $\theta_0=(0.4,0.6,0.7)$，迭代后得到模型参数 $\hat \theta=(0.4,0.537,0.643)$，说明EM算法与初始值的选择有关。

2 EM算法

2.1 算法描述

• 输入：观测变量数据 $X=(x^{(1)},x^{(2)},...x^{(m)})$，隐变量数据 $Z=(z^{(1)},z^{(2)},...z^{(m)})$。
• 输出：模型参数 $\theta$

观测数据 $X$又称不完全数据（incomplete-data），其概率分布是 $P(X \mid \theta)$
$X$和 $Z$连在一起称完全数据（complete data），其联合概率分布是 $P(Y,Z \mid \theta)$

• 优化目标：
  最大化观测数据 $X$关于参数 $\theta$的对数似然函数:
  $\begin{aligned} L(\theta) &= \log P(X \mid \theta) \\ &= \sum_X \log P(x \mid \theta) \\ &= \sum_X \log \sum_Z P(x,z \mid \theta ) \\ &= \sum_X \log \sum_Z P(z \mid \theta )P(x \mid z,\theta) \end{aligned}$即 $\begin{aligned} \hat \theta &= \arg \max_\theta L(\theta) \\ &= \arg \max_\theta\sum_X \log \sum_Z P(x,z \mid \theta ) \end{aligned}$

2.2 推导过程

$\begin{aligned} L(\theta) = \sum_X \log \sum_Z P(x,z \mid \theta ) &= \sum\limits_{i=1}^m \log\sum\limits_{z^{(i)}}P(x^{(i)},z^{(i)};\theta) \\ & = \sum\limits_{i=1}^m \log\sum\limits_{z^{(i)}}Q_i(z^{(i)})\frac{P(x^{(i)}， z^{(i)};\theta)}{Q_i(z^{(i)})} \\ & \geq \sum\limits_{i=1}^m \sum\limits_{z^{(i)}}Q_i(z^{(i)})\log\frac{P(x^{(i)}， z^{(i)};\theta)}{Q_i(z^{(i)})} \end{aligned}$令 $\frac{P(x^{(i)}， z^{(i)};\theta)}{Q_i(z^{(i)})} =c \;\;\; c为常数$由于 $Q_i(z^{(i)})$是一个分布，满足 $Q_i(z^{(i)})>0, \;\;\; \sum\limits_{z}Q_i(z^{(i)}) =1$则
$\begin{aligned} & \sum\limits_{z}P(x^{(i)}， z^{(i)};\theta) = c\sum\limits_{z}Q_i(z^{(i)}) \\ \Rightarrow & \;\;\; P(x^{(i)};\theta) =c \\ \Rightarrow & \;\;\; Q_i(z^{(i)}) = \frac{P(x^{(i)},z^{(i)};\theta)}{c}= \frac{P(x^{(i)},z^{(i)};\theta)}{P(x^{(i)};\theta)} = P( z^{(i)} \mid x^{(i)};\theta)) \end{aligned}$

通过极大化包含隐藏数据的对数似然的下界，从而极大化对数似然$ L(\theta)$。

此时，优化目标成为: $\arg \max \limits_{\theta} \sum\limits_{i=1}^m \sum\limits_{z^{(i)}}Q_i(z^{(i)})\log\frac{P(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}$去掉常数部分
$\arg \max \limits_{\theta} \sum\limits_{i=1}^m \sum\limits_{z^{(i)}}Q_i(z^{(i)})\log{P(x^{(i)},z^{(i)};\theta)}$
上式也就是EM算法的M步。

E步则是 $\log P(x^{(i)},z^{(i)};\theta)$基于条件概率分布 $Q_i(z^{(i)})$的期望。

2.3 算法步骤

EM通过迭代的方式逐步逼近 $L(\theta)$的极大值（反复构造新的下界）： $j$次迭代后 $\theta$的估计值为 $\theta_j$，希望新的估计值 $\theta_{j+1}$能使 $L(\theta)$的值增加，即 $L(\theta_{j+1}) > L(\theta_j)$。

具体步骤如下：

• 选泽参数初始值 $\theta_0$。
• E步：计算联合分布的条件概率期望，完全数据的对数似然函数 $\log P(X,Z \mid \theta)$在给定观测数据 $X$和当前参数 $\theta_j$下对未观测数据 $Z$的条件概率分布 $P(Z \mid X,\theta_j)$的期望：
  $\begin{aligned}\mathcal L(\theta,\theta_j) &= E_Z[\log P(X,Z \mid \theta) \mid X,\theta_j] \\ \\ &= \sum_X \sum_Z P(z \mid x,\theta_j) \log P(x,z \mid \theta) \\ &= \sum\limits_{i=1}^m\sum\limits_{z^{(i)}}Q_i(z^{(i)})\log{P(x^{(i)},z^{(i)};\theta)} \end{aligned}$其中， $Q_i(z^{(i)}) = P( z^{(i)} \mid x^{(i)},\theta^{j}))$
• M步：求第 $j+1$次迭代的参数估计值 $\theta_{j+1}$： $\theta_{j+1} = \arg \max_\theta \mathcal L(\theta,\theta_j)$
• 重复EM步，直至满足收敛条件，
  一般给出较小的正数 $\varepsilon_1,\varepsilon_2$，若满足
  $\mid\mid \theta_{j+1}-\theta_j \mid\mid \leq \varepsilon_1 或 \mid\mid \mathcal L(\theta_{j+1},\theta_{j}) - \mathcal L(\theta_{j+1},\theta_{j})\mid\mid \leq \varepsilon_2$则停止迭代，输出此时的 $\theta$值。

$\mathcal L$函数是EM算法的核心（有的教材中称其为 $Q$函数），每次迭代实际在求期望（E-步）及其极大化（M-步），实现参数 $\theta_j \rightarrow \theta_{j+1}$的更新，使似然函数增大或达到局部极值。

另一种理解

EM算法还可以看做用坐标下降法来最大化对数似然下界的过程，从而得到对数似然函数取得极大值时的参数估计，也就是极大似然估计的过程。
在迭代的过程中，每次固定一个变量，对另外的求极值，最后逐步逼近极值：

• E步：固定 $\theta$，优化 $\mathcal L$
• M步：固定 $\mathcal L$，优化 $\theta$
  交替将极值推向最大。

2.4 算法的收敛性

可以证明：
$\sum\limits_{i=1}^m \log P(x^{(i)};\theta_{j+1}) \geq \sum\limits_{i=1}^m \log P(x^{(i)};\theta_{j})$

EM算法可以收敛到局部极大值，但是却不能保证收敛到全局的极大值，因此它是局部最优的算法。当然，如果我们的优化目标 $L(\theta,\theta_{j})$是凸的，则EM算法可以保证收敛到全局最大值，这点和梯度下降法相同。

3 EM算法与非监督学习

• 监督学习中，训练数据 ${(x^{(1)},y^{(1)}),(x^{(2)},y^{(2)}),...,(x^{(m)},y^{(m)})}$中的每个样本点由输入和输出对组成，通过学习决策函数 $y=f(x)$（判别模型）或学习条件概率分布 $P(y \mid x)$（生成模型），实现数据的分类及回归。
• 非监督学习中，训练数据 ${(x^{(1)}, ),(x^{(2)},),...,(x^{(m)},)}$中每个样本点只有输入没有对应的输出，如聚类。

EM算法可以用于生成模型的非监督学习，生成模型由联合概率分布 $P(X,Y)$表示，认为非监督学习训练数据是联合概率分布产生的数据。

EM算法是一种框架，逼近统计模型参数的最大似然或最大后验估计。

K-Means可以看做一种Hard EM算法，

• E步：给定当前簇中心，每个对象都被指派到簇中心离该对象最近的簇（期望每个对象都属于最近的簇）。
• M步：给定指定簇，对于每个簇，算法调整其中心，使得指派到该簇的对象到该中心的距离之和最小化（指派到一个簇的对象的相似度最大化）。

K-Means是模糊聚类的一种特例（隐变量是存在分布的），在模糊或基于概率模型的聚类，EM算法从初始参数出发，并且迭代直到不能改善聚类（聚类收敛或改变充分小）。

• E-步：根据当前模糊聚类或概率簇的参数，把对象指派到簇中。
• M-步：发现新的聚类或参数，最小化模糊聚类的SSE或基于概率模型的期望似然。

4 EM算法应用于高斯混合模型学习

4.1 高斯混合模型

高斯混合模型（Gaussian misture model）指具有如下形式的概率分布模型：
$P(x \mid \theta) = \sum_{j=1}^{k} \alpha_j \phi(x \mid \theta_j)$
其中， $\alpha_j \geq 0$是系数， $\sum_{j=1}^{k} \alpha_j = 1$；
$\phi(x \mid \theta_j)$为高斯分布密度， $\theta_j = (\mu_j ,\theta^2_j)$
$\phi(x \mid \theta_j) = \frac{1}{\sqrt{2\pi}\sigma_j} \exp (-\frac{(x-\mu_j)^2}{2\sigma^2_j})$
称为第 $j$个分模型。

4.2 高斯混合模型参数估计

假设观测数据 ${x_1,x_2,...,x_m}$高斯混合模型 $P(x \mid \theta) = \sum_{j=1}^{k} \alpha_j \phi(x \mid \theta_j)$生成，其中模型参数 $\theta = (\alpha_1,\alpha_2,...,\alpha_k;\theta_1,\theta_2,...,\theta_k)$。

可以设想观测数据 ${x_1,x_2,...,x_m}$产生过程如下：

• 首先根据概率 $\alpha_j$选择第 $j$个高斯分布模型 $\phi(x \mid \theta_j)$；
• 然后根据第 $j$个模型的概率密度函数 $\phi(x \mid \theta_j)$生成观测数据 $x_i$。

执行上述过程 $m$次来生成观测数据。

反映观测数据 $x_i$来自第 $j$个模型的隐变量 $\gamma$是未知的，定义如下：
$\gamma_{ij} = \begin{cases} 1， & 第i个观测来自第j个分模型 \\ 0， & otherwise \end{cases} \;\; i=1,2,...,m;j=1,2,...,k$
$\gamma_{ij}$是0-1随机变量。

• E步：根据当前模型参数，计算模型 $i$对观测数据 $y_i$的响应度：
  $\hat \gamma_{ij} = \frac {\alpha_j \phi(x_i \mid \theta_j)}{\sum_{j=1}^{k}\alpha_j \phi(x_i \mid \theta_j)}$
• M步：计算新一轮迭代的模型参数：
  $\hat \mu_j = \frac{\sum_{i=1}^{m}\hat \gamma_{ij}y_i}{\sum_{i=1}^{m}\hat \gamma_{ij}}$

$\hat \theta_j^2 = \frac{\sum_{i=1}^{m}\hat \gamma_{ij}(x_i-\mu_j)^2}{\sum_{i=1}^{m}\hat \gamma_{ij}}$

$\hat \alpha_j = \frac{\sum_{i=1}^{m}\hat \gamma_{ij}}{m}$

EM算法还可以解释为F函数（F function）的极大-极大算法（maximization-maximization algorithm），基于这个解释有若干变形与推广，如广义期望极大算法（GEM，Generalized Expectation Maximization）

附：数学基础

（1）Jensen不等式

• 如果 $f$是凸函数， $X$是随机变量，则 $f(E[X]) \leqslant E[f(X)]$，特别地，如果 $f$是严格凸函数，那么等号成立当且仅当 $P(x=E[X])=1$时（ $X$是常量）；
• 如果 $f$是凹函数， $X$是随机变量，则 $f(E[X]) \geqslant E[f(X)]$，特别地，如果 $f$是严格凹函数，那么等号成立当且仅当 $P(x=E[X])=1$时（ $X$是常量）；
• 当 $f$是（严格）凹函数当且仅当 $-f$是（严格）凸函数。
  
  函数 $f$是凸函数， $X$是随机变量，有0.5的概率是 $a$，有0.5的概率是 $b$。 $X$的期望值 $E(x)$就是 $(a+b)/2$，可以看出 $E[f(X)]\geqslant f(E[X])$成立。

（2） 坐标下降法

坐标下降法是一种对多个参数分步迭代最终得到目标函数局部极大值（极小值）的方法。
如求目标函数 $f(\mu,\theta)$极值，可给定初值 $\theta_0$，然后分别计算：
$\mu_1= \arg \max_\mu f\left(\mu,\theta_0\right)$ $\theta_1=\arg \max_\theta f\left(\mu_1,\theta\right)$ $\mu_2=\arg \max_\mu f\left(\mu,\theta_1\right)$ $\theta_2=\arg \max_\theta f\left(\mu_2,\theta\right)$显然， $f\left( \mu_1,\theta_0 \right) \leqslant f\left( \mu_2,\theta_1 \right)... \leqslant f\left( \mu_{i+1},\theta_i \right)\leqslant f_{max}$

不断迭代后， $f\left( \mu_{i+1},\theta_i \right)$的值也逐渐逼近 $f_{max}$，此时 $\theta_i$和 $\mu_{i+1}$即为参数估计的值。