使用高斯混合模型(GMM)近似未知分布：EM算法的应用

该篇博客是对邱锡鹏老师《神经网络与深度学习》的学习笔记。在阅读本博文之前，建议读者先阅读上一篇博客EM算法。

高斯混合模型(Gaussian Mixture Model)

如果一个连续随机变量或连续随机向量的分布比较复杂，那么我们通常可以用高斯混合模型来估计其分布情况。

不失一般性，我们考虑一维的情况。假设样本 $x$ 是由$K$个高斯分布中的一个分布生成的，但是无法观测到具体由哪个分布生成。我们引入一个隐变量 z ∈ { 1 , ⋯   , K } z \in \{1,\cdots,K\} z∈{1,⋯,K}来表示样本 $x$ 来自于哪个高斯分布，$z$服从多项分布：
$$p(z=k;\pi )={\pi }_k,1\le k\le K$$

其中$\pi =[{\pi }_1,\cdots ,{\pi }_K]$为多项分布的参数，并满足${\pi }_k\ge 0,\forall k$，${\sum }_{k=1}^K{\pi }_k=1$，${\pi }_k$表示样本 $x$ 由第$k$个高斯分布生成的概率。给定$z=k$，条件分布$p(x|z=k)$为高斯分布：
$$\begin{array}{rl}p(x|z=k;{\mu }_k,{\sigma }_k)& =\mathcal{N}(x;{\mu }_k,{\sigma }_k)\\ & =\frac{1}{\sqrt{2\pi }{\sigma }_k}\exp \left(-\frac{(x-{\mu }_k{)}^2}{2{\sigma }_k^2}\right)\end{array}$$

其中${\mu }_k$和${\sigma }_k$表示第$k$个高斯分布的均值和方差。

从高斯混合模型GMM中生成一个样本 $x$ 的过程可以分为两步：
(1) 首先根据多项分布$p(z;\pi )$随机选取一个高斯分布；
(2) 假设选中第$k$个高斯分布（即$z=k$），再从高斯分布$\mathcal{N}(x;{\mu }_k,{\sigma }_k)$中选取一个样本$x$。

下图给出了高斯混合模型的图模型表示：

在GMM中，随机变量$x$的PDF为：
$$p(x)=\sum _{k=1}^K\mathcal{N}(x;{\mu }_k.{\sigma }_k)$$

参数估计

给定$N$个由高斯混合模型生成的训练样本$x^{(1)},x^{(2)},\cdots ,x^{(N)}$，希望能学习其中的参数${\pi }_k,{\mu }_k,{\sigma }_k,1\le k\le K$。由于我们无法观测样本$x^{(n)}$是从哪个高斯分布生成的，因此无法直接用最大似然来进行参数估计。

对每个样本$x^{(n)}$，其对数边际分布为
$$\log p(x^{(n)})=\log \sum _{k=1}^K{\pi }_k\cdot \mathcal{N}(x^{(n)};{\mu }_k,{\sigma }_k)$$

根据EM算法，参数估计分为两步进行迭代
(1) E步：先固定参数$\mu ,\sigma$，计算后验分布$p(z^{(n)}|x^{(n)})$，即
$$\begin{array}{rl}{\gamma }_{nk}& \doteq p(z^{(n)}=k|x^{(n)})\\ & =\frac{p(z^{(n)})p(x^{(n)}|z^{(n)})}{p(x^{(n)})}\\ & =\frac{{\pi }_k\cdot \mathcal{N}(x^{(n)};{\mu }_k,{\sigma }_k)}{{\sum }_{k=1}^K{\pi }_k\cdot \mathcal{N}(x^{(n)};{\mu }_k,{\sigma }_k)}\end{array}$$

其中${\gamma }_{nk}$定义了样本$x^{(n)}$属于第$k$个高斯样本的后验概率。

(2) M步：固定${\gamma }_{nk}$，令$q(z=k)={\gamma }_{nk}$，那么数据集$\mathcal{D}$的证据下界ELBO为：
$$\begin{array}{rl}ELBO(\gamma ,\mathcal{D};\pi ,\mu ,\sigma )& =\sum _{n=1}^N\sum _{k=1}^K{\gamma }_{nk}\log \frac{p(x^{(n)},z^{(n)}=k)}{{\gamma }_{nk}}\\ & =\sum _{n=1}^N\sum _{k=1}^K{\gamma }_{nk}\log \frac{p(z^{(n)}=k)\cdot p(x^{(n)}|z^{(n)}=k)}{{\gamma }_{nk}}\\ & =\sum _{n=1}^N\sum _{k=1}^K{\gamma }_{nk}\left(\log \mathcal{N}(x^{(n)};{\mu }_k,{\sigma }_k)+\log \frac{{\pi }_k}{{\gamma }_{nk}}\right)\end{array}$$

注意，上式中的${\gamma }_{nk}$的表达式中虽然含有参数${\pi }_k,{\mu }_k,{\sigma }_k$，但是这些参数都是第$t$次迭代估计的结果，即${\gamma }_{nk}$（后验），且$q(z=k)={\gamma }_{nk}$，在M步中被固定。而第M步要估计的第$t+1$次迭代的参数，因此，我们可以进一步把上式化简为：
$$ELBO(\gamma ,\mathcal{D};\pi ,\mu ,\sigma )=\sum _{n=1}^N\sum _{k=1}^K{\gamma }_{nk}\left(-\frac{(x-{\mu }_k{)}^2}{2{\sigma }_k^2}-\log {\sigma }_k+\log {\pi }_k\right)+Const$$

其中$Const$是和第$t+1$次M步迭代参数无关的常数。进一步将参数估计问题转化为优化问题：
$$\begin{array}{rl}& \underset{\pi ,\mu ,\sigma }{\max }ELBO(\gamma ,\mathcal{D};\pi ,\mu ,\sigma )\\ & s.t.\sum _{k=1}^K{\pi }_k=1\end{array}$$

利用拉格朗日乘数法求解上面的等式约束优化问题，分别求拉格朗日函数$ELBO(\gamma ,\mathcal{D};\pi ,\mu ,\sigma )+\lambda \left({\sum }_{k=1}^K{\pi }_k-1\right)$关于${\pi }_k,{\mu }_k,{\sigma }_k$的偏导数，令其等于0，可得
$$\begin{array}{rl}{\pi }_k& =\frac{N_k}{N}\\ {\mu }_k& =\frac{1}{N_k}\sum _{n=1}^N{\gamma }_{nk}{x}^{(n)}\\ {\sigma }_k^2& =\frac{1}{N_k}\sum _{n=1}^N{\gamma }_{nk}{\left(x^{(n)}-{\mu }_k\right)}^2\end{array}$$

其中
$$N_k=\sum _{n=1}^N{\gamma }_{nk}$$

参考

[1] 邱锡鹏，神经网络与深度学习，机械工业出版社，https://nndl.github.io/, 2020.