EM算法用于含有隐变量的概率模型参数的极大似然估计。这里首先对隐变量解释，举下面的例子

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

其中掷硬币 $A$ 的结果是未观测的，叫做隐变量，记做 $z$ 。将观测数据表示为 $Y=(Y_1,Y_2,\cdots,Y_n)^T$ ，未观测数据表示为 $Z=(Z_1,Z_2,\cdots,Z_n)^T$ ，则观测数据的似然函数为
$P(Y|\theta)=\sum_ZP(Y,Z|\theta)=\sum_ZP(Z|\theta)P(Y|Z,\theta)\tag1$
即
$P(Y|\theta)=\prod_{j=1}^n[\pi p^{y_j}(1-p)^{1-y_i}+(1-\pi)q^{y_j}(1-q)^{1-y_j}]\tag2$
对模型参数 $\theta=(\pi,p,q)$ 进行极大似然估计，即
$\hat\theta=\arg\max_\theta\log P(Y|\theta)\tag3$
因为掷硬币 $A$ 的结果未知，所以没有这个问题解析解，只能通过迭代的方式，逐步增加对数似然函数，找到一个解。EM算法解决的就是这样的一类问题。

接下来首先提出EM算法，然后对其进行解释。

EM算法

EM算法叫做Exception maximization算法，顾名思义，包含求期望（Exception ）和极大化（maximization）两个步骤。

输入：观测变量数据 $Y$ ，隐变量数据 $Z$ ，联合分布 $P(Y,Z|\theta)$ ，条件分布 $P(Z|Y,\theta)$ ；

输出：模型参数 $\theta$ 。

（1）选择参数的初值 $\theta^{(0)}$ ，开始迭代；

（2）E步：记 $\theta^{(i)}$ 为第 $i$ 次迭代参数 $\theta$ 的估计值，在第 $i+1$ 次迭代的E步，计算
$\begin{align} Q(\theta,\theta^{(i)})&=E_Z[\log P(Y,Z|\theta)|Y,\theta^{(i)}]\\ &=\sum_ZP(Z|Y,\theta^{(i)})\log P(Y,Z|\theta) \end{align}$
这里 $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\max_\theta Q(\theta,\theta^{(i)})$
（4）重复第(2)步和第(3)步，直到收敛。

第(2)步中 $Q(\theta,\theta^{(i)})$ 是EM算法的核心，称为Q函数。

Q函数定义：完全数据的对数似然函数 $\log P(Y,Z|\theta)$ 关于给定观测数据 $Y$ 和当前参数 $\theta^{(i)}$ 下对未观测数据 $Z$ 的条件概率分布 $P(Z|Y,\theta^{(i)})$ 的期望称为Q函数，即
$Q(\theta,\theta^{(i)})=E_Z[\log P(Y,Z|\theta)|Y,\theta^{(i)}]$
因为 $E_Z[g(z)]=\sum_Zp(z)g(z)$ ，所以
$\begin{align} Q(\theta,\theta^{(i)})&=E_Z[\log P(Y,Z|\theta)|Y,\theta^{(i)}]\\ &=\sum_ZP(Z|Y,\theta^{(i)})\log P(Y,Z|\theta) \end{align}\tag4$
其中，因为数据 $Y$ 未包含隐变量 $Z$ 的结果，所以称为不完全数据， $\log P(Y|\theta)$ 称为不完全数据的对数似然函数，而数据 $Y,Z$ 则称为完全数据， $\log P(Y,Z|\theta)$ 称为完全数据的对数似然函数。
下面关于EM算法作几点说明：

步骤（1）参数的初值可以任意选择，但需注意EM算法对初值是敏感的。
步骤（2）E步求 $Q(\theta,\theta^{(i)})$ 。Q函数式中 $Z$ 是未观测数据， $Y$ 是观测数据，注意， $Q(\theta,\theta^{(i)})$ 的第1个变元表示要极大化的参数，第2个变元表示参数的当前估计值。每次迭代实际在求Q函数及其极大。
步骤（3）M步求 $Q(\theta,\theta^{(i)})$ 的极大化，得到 $\theta^{(i+1)}$ ，完成一次迭代 $\theta^{(i)}\rightarrow\theta^{(i+1)}$ 。
步骤（4）给出停止迭代的条件，一般是对较小的正数 $\varepsilon_1$ ， $\varepsilon_2$ ，若满足

$||\theta^{(i+1)}-\theta^{(i)}||\lt\varepsilon_1\ \ 或\ \ ||Q(\theta^{(i+1)},\theta^{(i)})-Q(\theta^{(i)},\theta^{(i)})||\lt\varepsilon_2$

则停止迭代。

下面给出这种做法为什么可以对观测数据（不完全数据）进行极大似然估计。

EM算法的导出

但是如果通过第 $i$ 次迭代得到估计的参数 $\theta^{(i)}$ ，此时再找到一个参数 $\theta$ ，使得 $L(\theta)\gt L(\theta^{(i)})$ （和IIS算法思路有点类似），那么同样也可以起到极大似然估计的效果。为此，考虑两者的差
$L(\theta)-L(\theta^{(i)})=\log\bigg(\sum_ZP(Y|Z,\theta)P(z|\theta)\bigg)-\log P(Y|\theta^{(i)})$

利用Jensen不等式，过程略，得到其下界：
$L(\theta)-L(\theta^{(i)})\geq\sum_ZP(Z|Y,\theta^{(i)})\log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}\tag6$
于是得到对数似然函数的下界
$L(\theta)\geq L(\theta^{(i)})+\sum_ZP(Z|Y,\theta^{(i)})\log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}\tag7$
定义
$B(\theta,\theta^{(i)})\hat=L(\theta^{(i)})+\sum_ZP(Z|Y,\theta^{(i)})\log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}\tag8$
则
$L(\theta)\geq B(\theta,\theta^{(i)})$
并且
$\begin{align} B(\theta^{(i)},\theta^{(i)})&=L(\theta^{(i)})+\sum_ZP(Z|Y,\theta^{(i)})\log\frac{P(Y|Z,\theta^{(i)})P(Z|\theta^{(i)})}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}\\ &=L(\theta^{(i)})+\sum_ZP(Z|Y,\theta^{(i)})\log\frac{P(Y,Z|\theta^{(i)})}{P(Y,Z|\theta^{(i)})}\\ &=L(\theta^{(i)}) \end{align}$
即可以使 $B(\theta,\theta^{(i)})$ 可以增大的 $\theta$ ，也可以使 $L(\theta)$ 增大，所以极大似然估计变成了使 $B(\theta,\theta^{(i)})$ 极大化的问题，即
$\theta^{(i+1)}=\arg\max_\theta B(\theta,\theta^{(i)})$
所以
$\theta^{(i+1)}=\arg\max_\theta B(\theta,\theta^{(i)})=\arg\max_\theta\bigg(L(\theta^{(i)})+\sum_ZP(Z|Y,\theta^{(i)})\log\frac{P(Y|Z,\theta)P(Z|\theta)}{P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})}\bigg)$
省去不包含变量 $\theta$ 的常数项 $L(\theta^{(i)})$ ， $P(Z|Y,\theta^{(i)})P(Y|\theta^{(i)})$ ，得到
$\begin{align} \theta^{(i+1)}&=\arg\max_\theta\bigg(\sum_ZP(Z|Y,\theta^{(i)})\log P(Y|Z,\theta)P(Z|\theta)\bigg)\\ &=\arg\max_\theta\bigg(\sum_ZP(Z|Y,\theta^{(i)})\log P(Y,Z|\theta)\bigg)\\ &=\arg\max_\theta Q(\theta,\theta^{(i)}) \end{align}$
所以，在EM算法中最大化Q函数，就等同于最大对数似然函数的下界，从而进行极大似然估计。但是这种算法并不能保证找到全局最优值。

EM算法可以用于无监督学习，略。EM算法的收敛性证明略。

EM算法在高斯混合模型学习中的应用

高斯混合模型定义

高斯混合模型是指具有如下形式的概率分布模型：
$P(y|\theta)=\sum_{k=1}^K\alpha_k\phi(y|\theta_k)\tag9$
其中， $\alpha_k$ 是系数， $\alpha_k\geq0$ ， $\sum_{k=1}^K\alpha_k=1$ ； $\phi(y|\theta_k)$ 是高斯分布密度， $\theta_k=(\mu_k,\sigma_k^2)$ ，
$\phi(y|\theta_k)=\frac{1}{\sqrt{2\pi}\sigma_k}\exp\bigg(-\frac{(y-\mu_k)^2}{2\sigma_k^2}\bigg)$
称为第 $k$ 个模型。

高斯混合模型参数估计的EM算法

可以设想观测数据 $y_j$ ， $j=1,2,\cdots,N$ ，是这样产生的：首先，依概率 $\alpha_k$ 选择第 $k$ 个高斯分布分模型 $\phi(y|\theta_k)$ ；然后依第 $k$ 个分模型的概率分布 $\phi(y|\theta_k)$ 生成观测数据 $y_j$ ，这是观测数据 $y_j$ ， $j=1,2,\cdots,N$ ，是已知的；反应观测数据 $y_j$ 来自第 $k$ 个分模型的数据是未知的， $k=1,2,\cdots,K$ ，以隐变量 $\gamma_{jk}$ 表示，其定义如下：
$\gamma_{jk}= \begin{cases}1,\ \ &第j个观测来自第k个分模型\\ 0,\ \ & 否则 \end{cases}\\ j=1,2,\cdots,N;\ \ k=1,2,\cdots,K$
$\gamma_{jk}$ 是0,1随机变量。

有了观测数据 $y_j$ 及未观测数据 $\gamma_{jk}$ ，那么完全数据是
$(y_j,\gamma_{j1},\gamma_{j2},\cdots,\gamma_{jK}),\ \ j=1,2,\cdots,N$
于是可以写出完全数据的似然函数：
$\begin{align} P(y,\gamma|\theta)&=\prod_{j=1}^NP(y_j,\gamma_{j1},\gamma_{j2},\cdots,\gamma_{jK}|\theta)\\ &=\prod_{k=1}^K\prod_{j=1}^N[\alpha_k\phi(y_j|\theta_k)]^{\gamma_{jk}}\\ &=\prod_{k=1}^K\alpha_k^{n_k}\prod_{j=1}^N[\phi(y_j|\theta_k)]^{\gamma_{jk}}\\ &=\prod_{k=1}^K\alpha_k^{n_k}\prod_{j=1}^N\bigg[\frac{1}{\sqrt{2\pi}\sigma_k}\exp\bigg(-\frac{(y-\mu_k)^2}{2\sigma_k^2}\bigg)\bigg]^{\gamma_{jk}} \end{align}$
其中， $n_k=\sum_{j=1}^N\gamma_{jk}$ （依赖第 $k$ 个分类器的样本数）， $\sum_{k=1}^Kn_k=N$ 。

那么完全数据的对数似然函数为
$\log P(y,\gamma|\theta)=\sum_{k=1}^K\bigg\{n_k\log\alpha_k+\sum_{j=1}^N\gamma_{jk}\bigg[\log(\frac1{\sqrt{2\pi}})-\log\sigma_k-\frac1{2\sigma_k^2}(y_j-\mu_k^2)\bigg]\bigg\}$
确定Q函数
$\begin{align} Q(\theta,\theta^{(i)})&=E_{P(\gamma|y,\theta^{(i)})}[\log P(y,\gamma|\theta)]\\ &=E_{P(\gamma|y,\theta^{(i)})}\sum_{k=1}^K\bigg\{n_k\log\alpha_k+\sum_{j=1}^N\gamma_{jk}\bigg[\log(\frac1{\sqrt{2\pi}})-\log\sigma_k-\frac1{2\sigma_k^2}(y_j-\mu_k^2)\bigg]\bigg\}\\ &=\sum_{k=1}^K\bigg\{\sum_{j=1}^N(E\gamma_{jk})\log\alpha_k+\sum_{j=1}^N(E\gamma_{jk})\bigg[\log(\frac1{\sqrt{2\pi}})-\log\sigma_k-\frac1{2\sigma_k^2}(y_j-\mu_k^2)\bigg]\bigg\}\\ \end{align}\tag{10}$
其中
$\begin{align} E\gamma_{jk}&=\hat\gamma_{jk}=E(\gamma_{jk}|y,\theta)=P(\gamma_{jk}=1|y,\theta)\\ &=\frac{P(\gamma_{jk}=1,y_j|\theta)}{\sum_{k=1}^KP(\gamma_{jk}=1,y_j|\theta)}\\ &=\frac{P(y_j|\gamma_{jk}=1,\theta)P(\gamma_{jk}=1|\theta)}{\sum_{k=1}^KP(y_j|\gamma_{jk}=1,\theta)P(\gamma_{jk}=1|\theta)} \end{align}$

$\hat\gamma_{jk}$ 是在当前模型参数下第 $j$ 个观测数据来自第 $k$ 个分模型的概率，称为分模型 $k$ 对观测数据 $y_j$ 的响应度。

根据开头的描述， $P(y_j|\gamma_{jk}=1,\theta)=\phi(y_k|\theta_k)$ ， $P(\gamma_{jk}=1|\theta)=\alpha_k$ ，所以
$E\gamma_{jk}=\hat\gamma_{jk}=E(\gamma_{jk}|y,\theta)=\frac{\alpha_k\phi(y_j|\theta_k)}{\sum_{k=1}^K\alpha_k\phi(y_j|\theta_k)}\tag{11}$
将 $\hat\gamma_{jk}=E\gamma_{jk}$ 及 $n_k=\sum_{j=1}^NE\gamma_{jk}$ （此处的 $n_k$ 实际为 $E_{P(\gamma_{jk}|y,\theta)}n_k$ ， $n_k=\sum_{j=1}^N\gamma_{jk}$ ）代入公式(10)得到
$Q(\theta,\theta^{(i)})=\sum_{k=1}^K\bigg\{n_k\log\alpha_k+\sum_{j=1}^N\hat\gamma_{jk}\bigg[\log(\frac1{\sqrt{2\pi}})-\log\sigma_k-\frac1{2\sigma_k^2}(y_j-\mu_k^2)\bigg]\bigg\}\tag{12}$
接下来需要求Q函数的极大化（极大化观测数据对数似然函数的下界），即
$\theta^{(i+1)}=\arg\max_\theta Q(\theta,\theta^{(i)})$
其中 $\theta=(\alpha_k,\mu_k,\sigma_k^2)$ ， $k=1,2,\cdots,K$ ，将公式(12)对 $\mu_k$ 和 $\sigma_k^2$ 求偏导等于0，得到第 $i+1$ 次的更新值 $\hat\mu_k$ 和 $\hat\sigma_k^2$ ，同样将公式(12)对 $\alpha_k$ 求偏导等于0，加上条件 $\sum_{k=1}^K\alpha_k=1$ ，求得更新值 $\hat\alpha_k$ 。计算更新值时 $\hat\gamma_{jk}$ ， $n_k$ 用到的参数是第 $i$ 次更新得到的值，所以可以通过迭代的方式不断更新参数直到收敛。求得的 $\hat\mu_k$ ， $\hat\sigma_k^2$ 和 $\hat\alpha_k$ 为
$\hat\mu_k=\frac{\sum_{j=1}^N}{\sum_{j=1}^N\hat\gamma_{jk}},\ \ k=1,2,\cdots,K\tag{13}$

$\hat\sigma_k^2=\frac{\sum_{j=1}^N\hat\gamma_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N\hat\gamma_{jk}},\ \ k=1,2,\cdots,K\tag{14}$

$\hat\alpha_k=\frac{n_k}N=\frac{\sum_{j=1}^N\hat\gamma_{jk}}N,\ \ k=1,2,\cdots,K\tag{15}$

高斯混合模型参数估计的EM算法

输入：观测数据 $y_1,y_2,\cdots,y_N$ ，高斯混合模型；

输出：高斯混合模型参数。

（1）取参数的初始值开始迭代（初值敏感）

（2）E步：依据当前模型参数，计算分模型 $k$ 对观测数据 $y_j$ 的响应度
$\hat\gamma_{jk}=\frac{\alpha_k\phi(y_j|\theta_k)}{\sum_{k=1}^K\alpha_k\phi(y_j|\theta_k)},\ \ j=1,2,\cdots,N;\ \ k=1,2,\cdots,K$
（3）M步：计算新一轮迭代的模型参数
$\hat\mu_k=\frac{\sum_{j=1}^N}{\sum_{j=1}^N\hat\gamma_{jk}},\ \ k=1,2,\cdots,K$

$\hat\sigma_k^2=\frac{\sum_{j=1}^N\hat\gamma_{jk}(y_j-\mu_k)^2}{\sum_{j=1}^N\hat\gamma_{jk}},\ \ k=1,2,\cdots,K$

$\hat\alpha_k=\frac{n_k}N=\frac{\sum_{j=1}^N\hat\gamma_{jk}}N,\ \ k=1,2,\cdots,K$

（4）重复第(2)步和第(3)步，直到收敛。

统计机器学习-EM算法（期望极大算法）

统计机器学习-EM算法（期望极大算法）

EM算法

EM算法的导出

EM算法在高斯混合模型学习中的应用

高斯混合模型定义

高斯混合模型参数估计的EM算法

高斯混合模型参数估计的EM算法