EM求解的高斯混合模型——Q函数的极大似然估计（八）

先导：EM求解的混合密度模型——Q函数

$$p(\mathbf{x}\mid {\mathbf{\theta }}_k)\to N(\mathbf{x}\mid {\mathbf{\mu }}_{\mathbf{k}},{\mathbf{\Sigma }}_k)$$

由上述推导即可获得高斯混合模型的 EM 算法：在每步迭代中，先根据当前参数来计算每个样本属于每个高斯成分的后验概率${\gamma }_{ji}$（E步），再更新模型参数$\{({\pi }_i,{\mu }_i,{\Sigma }_i)\mid 1\le i\le k\}$（M步）。

---

E步：
$${\gamma }_k^{(t)}({\mathbf{x}}_i)=\frac{{\pi }_k^{(t)}N({\mathbf{x}}_i\mid {\mathbf{\mu }}_k^{(t)},{\Sigma }_k^{(t)})}{\sum _{k=1}^K{\pi }_k^{(t)}N({\mathbf{x}}_i\mid {\mathbf{\mu }}_k^{(t)},{\Sigma }_k^{(t)})}$$

---

可以发现，根据离散随机变量期望的定义，替换后的式子事实上是在对原先的对数似然函数计算期望（Q函数）：

$$Q(\theta \mid X,{\theta }^{(t)})=\sum _{j=1}^n\sum _{k=1}^K\gamma (z_{jk})\left(\ln {\pi }_k+\ln N_k({\mathbf{x}}_j\mid {\mathbf{\mu }}_k,{\Sigma }_k)\right)$$

其中，${\theta }^{(t)}$是当前的模型参数。根据式 (10)，上式中的$\gamma (z_{jk})$可以根据当前的模型参数计算出来，不涉及隐变量。

由于高斯混合模型是概率模型，无法肯定地指出某个样本属于某个聚类，因此，将隐变量$z_{ik}$替换为上文所介绍的责任$\gamma (z_{ik})=P(z_{ik}=1)=P(z_k=1|{\mathbf{x}}_i)$，即样本${\mathbf{x}}_i$属于第$k$个聚类的概率，这相当于一个软版本的二元隐变量。

---

其代入式 (18)，回顾 θ = { π k , μ k , Σ k } i = 1 K \theta = \{ \pi_k, {\bm \mu}_k, {\bm \varSigma}_k \}_{i=1}^{K} θ={πk​,μk​,Σk​}i=1K​，得到

Q ( { π k , μ k , Σ k } i = 1 K , { π k ( t ) , μ k ( t ) , Σ k ( t ) } i = 1 K ) = ∑ i = 1 n ∑ k = 1 K γ k ( t ) ( x i ) ln ⁡ ( π k N k ( x j ∣ μ k , Σ k ) ) (20) Q(\{ \pi_k, {\bm \mu}_k, {\bm \varSigma}_k \}_{i=1}^{K}, \{ \pi_k^{(t)}, {\bm \mu}_k^{(t)}, {\bm \varSigma}_k^{(t)} \}_{i=1}^{K}) \\= \sum\limits_{i=1}^{n} \sum\limits_{k=1}^{K} \gamma_k^{(t)}({\bm x}_i) \ln (\pi_k N_k(\bm{x}_j \mid {\bm \mu}_k, {\bm \varSigma}_k)) \tag{20} Q({πk​,μk​,Σk​}i=1K​,{πk(t)​,μk(t)​,Σk(t)​}i=1K​)=i=1∑n​k=1∑K​γk(t)​(xi​)ln(πk​Nk​(xj​∣μk​,Σk​))(20)

有了式 (19) 的期望对数似然函数，可以对其最大化，从而得到对模型参数的当前最优估计。

为了获得观测数据的聚类隶属度，需要找到一组参数 { π k , μ k , Σ k } i = 1 K \{ \pi_k, {\bm \mu}_k, {\bm \varSigma}_k \}_{i=1}^{K} {πk​,μk​,Σk​}i=1K​的估计。

---

M 步：

{ π k ( t + 1 ) , μ k ( t + 1 ) , Σ k ( t + 1 ) } i = 1 K = arg ⁡ max ⁡ { π k , μ k , Σ k } i = 1 K Q ( { π k , μ k , Σ k } i = 1 K , { π k ( t ) , μ k ( t ) , Σ k ( t ) } i = 1 K ) (21) \{ \pi_k^{(t+1)}, {\bm \mu}_k^{(t+1)}, {\bm \varSigma}_k^{(t+1)} \}_{i=1}^{K} \\= \arg \max_{\{ \pi_k, {\bm \mu}_k, {\bm \varSigma}_k \}_{i=1}^{K}} Q(\{ \pi_k, {\bm \mu}_k, {\bm \varSigma}_k \}_{i=1}^{K}, \{ \pi_k^{(t)}, {\bm \mu}_k^{(t)}, {\bm \varSigma}_k^{(t)} \}_{i=1}^{K}) \tag{21} {πk(t+1)​,μk(t+1)​,Σk(t+1)​}i=1K​=arg{πk​,μk​,Σk​}i=1K​max​Q({πk​,μk​,Σk​}i=1K​,{πk(t)​,μk(t)​,Σk(t)​}i=1K​)(21)

这个最大化是比较直接的。通过对${\mathbf{\mu }}_k$和${\Sigma }_k^{-1}$微分，有：

$$\frac{\partial Q(\theta \mid X,{\theta }^{(t)})}{\partial {\mathbf{\mu }}_k}=\sum _{i=1}^n{\gamma }_k^{(t)}({\mathbf{x}}_i){\Sigma }_k^{-1}({\mathbf{x}}_i-{\mathbf{\mu }}_k)=0$$

$$\frac{\partial Q(\theta \mid X,{\theta }^{(t)})}{\partial {\Sigma }_k^{-1}}=-\frac{1}{2}\sum _{i=1}^n{\gamma }_k^{(t)}({\mathbf{x}}_i)({\Sigma }_k-({\mathbf{x}}_i-{\mathbf{\mu }}_k)({\mathbf{x}}_i-{\mathbf{\mu }}_k{)}^{\top })=0\text{\hspace{1em}(22)}$$

其中$k=1,\cdots ,K$，解出方程组：

$${\mathbf{\mu }}_k^{(t+1)}=\frac{\sum _{i=1}^n{\gamma }_k^{(t)}({\mathbf{x}}_i){\mathbf{x}}_i}{\sum _{i=1}^n{\gamma }_k^{(t)}({\mathbf{x}}_i)}$$

$${\Sigma }_k^{(t+1)}=\frac{\sum _{i=1}^n{\gamma }_k^{(t)}({\mathbf{x}}_i)({\mathbf{x}}_i-{\hat{\mu }}_k^{(t+1)})({\mathbf{x}}_i-{\hat{\mu }}_k^{(t+1)}{)}^{\top }}{\sum _{i=1}^n{\gamma }_k^{(t)}({\mathbf{x}}_i)}\text{\hspace{1em}(23)}$$

其中$k=1,\cdots ,K$，可以证明这个驻点确实是一个极大值点。

至于${\pi }_k$的最大化，为约束$\sum _{k=1}^K{\pi }_k=1$，引入拉格朗日乘子$\lambda$从而寻找$Q(\theta \mid X,{\theta }^{(t)})+\lambda (\sum _{k=1}^K{\pi }_k-1)$的一个驻点。对${\pi }_k$微分得到：

$$\frac{\partial Q(\theta \mid X,{\theta }^{(t)})}{\partial {\pi }_k}=\sum _{i=1}^n\frac{{\gamma }_k^{(t)}({\mathbf{x}}_i)}{{\pi }_k}+\lambda =0\text{\hspace{1em}(24)}$$

其中$k=1,\cdots ,K$，把这 K 个方程的两边都乘以${\pi }_k$（假设所有的${\pi }_k$都不为零），然后把它们相加得到：

$$\sum _{i=1}^n\sum _{k=1}^K{\gamma }_k^{(t)}({\mathbf{x}}_i)+\lambda \sum _{k=1}^K{\pi }_k=0$$

由于$\sum _{k=1}^K{\gamma }_k^{(t)}({\mathbf{x}}_i)$和$\sum _{k=1}^K{\pi }_k$都等于 1，有
$$\lambda =-n\text{\hspace{1em}(25)}$$

把$\lambda =-n$代入到式 (24)，解方程得到

$${\pi }_k^{(t+1)}=\frac{1}{n}\sum _{i=1}^n{\gamma }_k^{(t)}({\mathbf{x}}_i)\text{\hspace{1em}(26)}$$

其中$k=1,\cdots ,K$。这个过程一直重复，直到对数似然$L(\theta )$没有显著变化。此时，如果需要，可以进行聚类分配，方法是将点${\mathbf{x}}_i$分配给拥有最大组成员隶属度的组$k$。注意，可以进行其他的硬聚类分配。例如，该算法可能拒绝为点分配任何组，除非最大的组隶属度超过预先指定的阈值。

整个 EM 过程总结如下。

---

算法 使用EM和高斯混合模型聚类

1. 初始化$K$，$\tau >0$， { π k ( 0 ) , μ k ( 0 ) , Σ k ( 0 ) } k = 1 K \{\pi_k^{(0)}, {\bm \mu}_k^{(0)}, {\bm \varSigma}_k^{(0)}\}_{k=1}^K {πk(0)​,μk(0)​,Σk(0)​}k=1K​

2. repeat

3. E步：更新组成员
   $${\gamma }_k^{(t)}({\mathbf{x}}_i)=\frac{{\pi }_k^{(t)}N({\mathbf{x}}_i\mid {\mathbf{\mu }}_k^{(t)},{\Sigma }_k^{(t)})}{\sum _{k=1}^K{\pi }_k^{(t)}N({\mathbf{x}}_i\mid {\mathbf{\mu }}_k^{(t)},{\Sigma }_k^{(t)})}$$

4. M步：重新估计模型参数
   $${\mathbf{\mu }}_k^{(t+1)}=\frac{\sum _{i=1}^n{\gamma }_k^{(t)}({\mathbf{x}}_i){\mathbf{x}}_i}{\sum _{i=1}^n{\gamma }_k^{(t)}({\mathbf{x}}_i)}$$$${\Sigma }_k^{(t+1)}=\frac{\sum _{i=1}^n{\gamma }_k^{(t)}({\mathbf{x}}_i)({\mathbf{x}}_i-{\hat{\mu }}_k^{(t+1)})({\mathbf{x}}_i-{\hat{\mu }}_k^{(t+1)}{)}^{\top }}{\sum _{i=1}^n{\gamma }_k^{(t)}({\mathbf{x}}_i)}$$$${\pi }_k^{(t+1)}=\frac{1}{n}\sum _{i=1}^n{\gamma }_k^{(t)}({\mathbf{x}}_i)$$

5. 计算对数似然：
   L ( { π k ( t + 1 ) , μ k ( t + 1 ) , Σ k ( t + 1 ) } k = 1 K ) = ∑ i = 1 n ln ⁡ ( ∑ k = 1 K π k ( t + 1 ) N ( x i ∣ μ k ( t + 1 ) , Σ k ( t + 1 ) ) ) L(\{\pi_k^{(t+1)}, {\bm \mu}_k^{(t+1)}, {\bm \varSigma}_k^{(t+1)}\}_{k=1}^K) = \sum\limits_{i=1}^n \ln \left( \sum\limits_{k=1}^K \pi_k^{(t+1)} N({\bm x}_i \mid {\bm \mu}_k^{(t+1)}, {\bm \varSigma}_k^{(t+1)}) \right) L({πk(t+1)​,μk(t+1)​,Σk(t+1)​}k=1K​)=i=1∑n​ln(k=1∑K​πk(t+1)​N(xi​∣μk(t+1)​,Σk(t+1)​))

6. until ∣ L ( { π k ( t + 1 ) , μ k ( t + 1 ) , Σ k ( t + 1 ) } k = 1 K ) − L ( { π k ( t ) , μ k ( t ) , Σ k ( t ) } k = 1 K ) ∣ < τ |L(\{\pi_k^{(t+1)}, {\bm \mu}_k^{(t+1)}, {\bm \varSigma}_k^{(t+1)}\}_{k=1}^K) - L(\{\pi_k^{(t)}, {\bm \mu}_k^{(t)}, {\bm \varSigma}_k^{(t)}\}_{k=1}^K)| < \tau ∣L({πk(t+1)​,μk(t+1)​,Σk(t+1)​}k=1K​)−L({πk(t)​,μk(t)​,Σk(t)​}k=1K​)∣<τ

---

总结

• 均值${\mu }_i$：通过样本加权平均来估计。
• 协方差矩阵${\Sigma }_i$：通过样本加权中心矩来估计。
• 混合系数${\pi }_i$：通过样本属于该成分的平均后验概率来确定。