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

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

高斯混合模型(Gaussian Mixture Model)

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

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

其中 π = [ π 1 , ⋯ , π K ] \pi = [\pi_1, \cdots, \pi_K] π=[π1,⋯,πK]为多项分布的参数，并满足 π k ≥ 0 , ∀ k \pi_k \geq 0, \forall k πk≥0,∀k， ∑ k = 1 K π k = 1 \sum_{k=1}^K \pi_k = 1 ∑k=1Kπk=1， π k \pi_k πk表示样本 x x x 由第 k k k个高斯分布生成的概率。给定 z = k z=k z=k，条件分布 p ( x ∣ z = k ) p(x|z=k) p(x∣z=k)为高斯分布：
p ( x ∣ z = k ; μ k , σ k ) = N ( x ; μ k , σ k ) = 1 2 π σ k exp ⁡ ( − ( x − μ k ) 2 2 σ k 2 ) \begin{aligned} 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^2_k} \right ) \end{aligned} p(x∣z=k;μk,σk)=N(x;μk,σk)=2π σk1exp(−2σk2(x−μk)2)

其中 μ k \mu_k μk和 σ k \sigma_k σk表示第 k k k个高斯分布的均值和方差。

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

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

在GMM中，随机变量 x x x的PDF为：
p ( x ) = ∑ k = 1 K N ( x ; μ k . σ k ) p(x) = \sum_{k=1}^K \mathcal{N} (x; \mu_k. \sigma_k) p(x)=k=1∑KN(x;μk.σk)

参数估计

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

对每个样本 x ( n ) x^{(n)} x(n)，其对数边际分布为
log ⁡ p ( x ( n ) ) = log ⁡ ∑ k = 1 K π k ⋅ N ( x ( n ) ; μ k , σ k ) \log p( x^{(n)}) = \log \sum_{k=1}^K \pi_k \cdot \mathcal{N}(x^{(n)}; \mu_k, \sigma_k) logp(x(n))=logk=1∑Kπk⋅N(x(n);μk,σk)

根据EM算法，参数估计分为两步进行迭代
(1) E步：先固定参数 μ , σ \mu, \sigma μ,σ，计算后验分布 p ( z ( n ) ∣ x ( n ) ) p(z^{(n)}|x^{(n)}) p(z(n)∣x(n))，即
γ n k ≐ p ( z ( n ) = k ∣ x ( n ) ) = p ( z ( n ) ) p ( x ( n ) ∣ z ( n ) ) p ( x ( n ) ) = π k ⋅ N ( x ( n ) ; μ k , σ k ) ∑ k = 1 K π k ⋅ N ( x ( n ) ; μ k , σ k ) \begin{aligned} \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{aligned} γnk≐p(z(n)=k∣x(n))=p(x(n))p(z(n))p(x(n)∣z(n))=∑k=1Kπk⋅N(x(n);μk,σk)πk⋅N(x(n);μk,σk)

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

(2) M步：固定 γ n k \gamma_{nk} γnk，令 q ( z = k ) = γ n k q(z=k)=\gamma_{nk} q(z=k)=γnk，那么数据集 D \mathcal D D的证据下界ELBO为：
E L B O ( γ , D ; π , μ , σ ) = ∑ n = 1 N ∑ k = 1 K γ n k log ⁡ p ( x ( n ) , z ( n ) = k ) γ n k = ∑ n = 1 N ∑ k = 1 K γ n k log ⁡ p ( z ( n ) = k ) ⋅ p ( x ( n ) ∣ z ( n ) = k ) γ n k = ∑ n = 1 N ∑ k = 1 K γ n k ( log ⁡ N ( x ( n ) ; μ k , σ k ) + log ⁡ π k γ n k ) \begin{aligned} 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{aligned} ELBO(γ,D;π,μ,σ)=n=1∑Nk=1∑Kγnklogγnkp(x(n),z(n)=k)=n=1∑Nk=1∑Kγnklogγnkp(z(n)=k)⋅p(x(n)∣z(n)=k)=n=1∑Nk=1∑Kγnk(logN(x(n);μk,σk)+logγnkπk)

注意，上式中的 γ n k \gamma_{nk} γnk的表达式中虽然含有参数 π k , μ k , σ k \pi_k, \mu_k, \sigma_k πk,μk,σk，但是这些参数都是第 t t t次迭代估计的结果，即 γ n k \gamma_{nk} γnk（后验），且 q ( z = k ) = γ n k q(z=k)=\gamma_{nk} q(z=k)=γnk，在M步中被固定。而第M步要估计的第 t + 1 t+1 t+1次迭代的参数，因此，我们可以进一步把上式化简为：
E L B O ( γ , D ; π , μ , σ ) = ∑ n = 1 N ∑ k = 1 K γ n k ( − ( x − μ k ) 2 2 σ k 2 − log ⁡ σ k + log ⁡ π k ) + C o n s t 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^2_k} - \log \sigma_k + \log \pi _k \right) + Const ELBO(γ,D;π,μ,σ)=n=1∑Nk=1∑Kγnk(−2σk2(x−μk)2−logσk+logπk)+Const

其中 C o n s t Const Const是和第 t + 1 t+1 t+1次M步迭代参数无关的常数。进一步将参数估计问题转化为优化问题：
max ⁡ π , μ , σ E L B O ( γ , D ; π , μ , σ ) s . t . ∑ k = 1 K π k = 1 \begin{aligned} & \max_{\pi, \mu, \sigma} ELBO(\gamma, \mathcal D; \pi, \mu, \sigma) \\ & s.t. \ \ \sum_{k=1}^K \pi_k = 1 \end{aligned} π,μ,σmaxELBO(γ,D;π,μ,σ)s.t. k=1∑Kπk=1

利用拉格朗日乘数法求解上面的等式约束优化问题，分别求拉格朗日函数 E L B O ( γ , D ; π , μ , σ ) + λ ( ∑ k = 1 K π k − 1 ) ELBO(\gamma, \mathcal D; \pi, \mu, \sigma) + \lambda \left ( \sum_{k=1}^K \pi_k -1 \right) ELBO(γ,D;π,μ,σ)+λ(∑k=1Kπk−1)关于 π k , μ k , σ k \pi_k, \mu_k, \sigma_k πk,μk,σk的偏导数，令其等于0，可得
π k = N k N μ k = 1 N k ∑ n = 1 N γ n k x ( n ) σ k 2 = 1 N k ∑ n = 1 N γ n k ( x ( n ) − μ k ) 2 \begin{aligned} \pi_k &= \frac{N_k}{N} \\ \mu_k &= \frac{1}{N_k} \sum_{n=1}^N \gamma_{nk} x^{(n)} \\ \sigma^2_k &= \frac{1}{N_k} \sum_{n=1}^N \gamma_{nk} \left ( x^{(n)} - \mu_k \right )^2 \\ \end{aligned} πkμkσk2=NNk=Nk1n=1∑Nγnkx(n)=Nk1n=1∑Nγnk(x(n)−μk)2

其中
N k = ∑ n = 1 N γ n k N_k = \sum_{n=1}^N \gamma_{nk} Nk=n=1∑Nγnk

参考

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