Expectation maximization algorithm

回顾一下GMM以及HMM等模型的求解方法，即著名的Expectation Maximization (EM) 算法。参考李航老师的《统计学习方法》，EM算法是一种迭代算法，用于含有隐含变量的概率模型参数的极大似然估计（MLE），或极大后验概率估计 (MPE).

EM算法

定义$Y$为观测随机变量，$Z$为隐含随机变量，则$Y$和$Z$一起构成完全数据 (complete-data)。假设待估计的概率模型参数为$\theta$，则观测数据的概率分布为$P(Y|\theta)$, 即为其似然函数，对应的对数似然函数为$L(\theta) = \log{P(Y|\theta)}$. 设$Y$和$Z$Z的联合概率分布为$P(Y,Z|\theta)$，那么完全数据的对数似然函数为$\log{P(Y,Z|\theta)}$.

EM算法的目的就是求解参数$\theta$，极大化观测量的对数似然函数$L(\theta)$。因为包含隐含量，所以采用迭代的方法，分为E (期望) 和M (极大化)两步，求解算法如下。

1. 初始化参数$\theta^{(0)}$;

2. E step: 记$\theta^{(i)}$为第$i$次迭代参数$\theta$的估计值，则第$i+1$次迭代，计算如下期望.

   $$\begin{align} Q(\theta,\theta^{(i)}) &= E_{z}[\log{P(Y,Z|\theta)}|Y,\theta^{(i)}] \notag \\ &=\sum_{z}{\log{P(Y,Z|\theta)P(Z|Y,\theta^{(i)})}}, \end{align}$$

   其中$P(Z|Y,\theta^{(i)})$是在给定观测数据$Y$和当前估计的参数$\theta^{(i)}$下隐含变量$Z$的条件概率分布；$Q(\theta,\theta^{(i)})$定义为完全数据的对数似然函数在给定观测数据和当前参数下对隐含数据$Z$的条件概率分布$P(Z|Y,\theta^{(i)})$的期望，即$\log{P(Y|\theta)}$。

3. M step: 求使$Q(\theta, \theta^{(i)})$极大化的参数$\theta$，以确定第$i+1$次迭代的参数的估计值$\theta^{(i+1)}$, 即计算

   $$\begin{align} \theta^{(i+1)} = \mathrm{arg}\max\limits_{\theta}{Q(\theta,\theta^{(i)})} \end{align}$$

4. 重复E和M两步，直到收敛.

李航老师也强调了一点：EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。 我在这个notebook里给了一个样例，即“三硬币模型”。

高斯混合模型的求解

EM算法最经典的应用就是高斯混合模型的参数估计，该模型是语音信号处理的基础模型之一，其定义如下，

$$\begin{align} P(\mathbf{y}|\mathbf{\theta}) = \sum_{k=1}^{K}{c_k\phi(\mathbf{y}|\mathbf{\theta_k)}} \end{align}$$

其中，$c_k$是系数，$c_k\geq0$, $\sum_{k=1}^{K}{c_k=1}$; $\phi(\mathbf{y}|\mathbf{\theta_k})$是高斯分布，$\mathbf{\theta_k} = (\mathbf{\mu_k},\Sigma_k^2)$,

$$\begin{equation} \phi(\mathbf{y}|\mathbf{\theta_k}) = \frac{1}{(2\pi)^{M/2}|\Sigma_k|}{-\exp{(\mathbf{y}-\mathbf{\mu_k})^{T}{\Sigma_{k}}^{-1}(\mathbf{y}-\mathbf{\mu_k})}}. \end{equation}$$

在GMM模型中有，

• 观测变量：$\mathbf{Y}$，
• 隐含变量：$\gamma_k\in{0,1}, k=1,~2,~,\cdots,~K$，表示$\mathbf{y}$是否来自第$k$个高斯分量
• 模型参数：$\mathbf{\theta} = {c_k, \mathbf{\theta_k}}, k=1,~2,~,\cdots,~K$.

详细的推导见《统计学习方法》，这里给出GMM模型的E和M步骤.
设观测数据$\mathbf{Y} = \{\mathbf{y1},~\mathbf{y2},~\cdots,~\mathbf{y_N}\}$

• E step

  $$\begin{align} \hat{\gamma_{jk}} = \frac{c_k\phi(\mathbf{y_j}|\mathbf{\theta_k})}{\sum_{k=1}^{K}{c_k\phi(\mathbf{y_j}|\mathbf{\theta_k})}} \end{align}$$

  其中$j=1,~2,~\cdots,~,N, k=1,~2,~\cdots,~K$.

• M step

  $$\begin{align} \hat{\mu}_{km} &= \sum_{j=1}^{N}{\hat{\gamma}_{jk}y_{jm}}/{\sum^{N}_{j=1}{\hat{\gamma}_{jk}}} \notag \\ \hat{\sigma}^{2}_{km} &= \sum_{j=1}^{N}{\hat{\gamma}_{jk}(y_{jm}-\mu_{km})^2}/{\sum_{j=1}^{N}{\hat{\gamma}_{jk}}} \notag \\ \hat{c}_{k} &= {\sum_{j=1}^{N}{\hat{\gamma}_{jk}}} / N \end{align}$$

  其中$m=1,~2,~\cdots,~,M, k=1,~2,~\cdots,~K$.

最后，给出一个样例，简单的二维GMM，如下图，notebook见这里.

可以看出，EM的估计效果还是不错的，并且提供初始化参数值的结果 (左下) 比随机初始化 (右下) 的结果要好。

Reference

[1] 李航，统计学系方法，2012 清华大学出版社
[2] Generating random variables from a mixture of Normal distributions