基于高斯混合模型(GMM)的期望最大化(EM)算法参数估计过程
题目描述
我们考虑一个数据集,其中的样本点可能来自多个不同的、但未直接观测到(即“隐变量”)的子群体(或“成分”)。每个子群体都服从一个高斯(正态)分布,但整体的数据分布是这些高斯分布的加权混合。高斯混合模型(Gaussian Mixture Model, GMM)就是用来描述这种混合分布的概率模型。我们的任务是:给定一组观测数据,如何估计出GMM中各个高斯成分的参数(均值、协方差)以及它们的混合权重?期望最大化(Expectation-Maximization, EM)算法为解决此类含有隐变量的概率模型的极大似然估计提供了一种优雅的迭代框架。
解题过程
我们以最大似然估计为目标,但模型中存在“隐变量”——即每个数据点具体来自于哪一个高斯成分。EM算法通过迭代执行两个步骤来绕过这个难题:E步(Expectation)基于当前参数估计隐变量的后验分布;M步(Maximization)基于E步计算出的“完全数据”的期望对数似然,更新模型参数以最大化这个期望。
步骤1:高斯混合模型的数学定义
假设我们有观测数据集 \(\mathbf{X} = \{ \mathbf{x}_1, \mathbf{x}_2, ..., \mathbf{x}_N \}\),其中 \(\mathbf{x}_i \in \mathbb{R}^D\)。GMM假设数据由 \(K\) 个高斯分布混合而成。模型参数为 \(\Theta = \{ \pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k \}_{k=1}^K\),其中:
- \(\pi_k\) 是第 \(k\) 个高斯成分的混合权重,满足 \(0 \leq \pi_k \leq 1\) 且 \(\sum_{k=1}^K \pi_k = 1\)。
- \(\boldsymbol{\mu}_k\) 是第 \(k\) 个高斯分布的均值向量。
- \(\boldsymbol{\Sigma}_k\) 是第 \(k\) 个高斯分布的协方差矩阵。
整个模型的概率密度函数是各个高斯分布密度的加权和:
\[p(\mathbf{x} | \Theta) = \sum_{k=1}^{K} \pi_k \mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \]
其中 \(\mathcal{N}(\mathbf{x} | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)\) 表示多元高斯分布的概率密度。
步骤2:引入隐变量与完全数据对数似然
为每个数据点 \(\mathbf{x}_i\) 引入一个 \(K\) 维的二元隐变量 \(\mathbf{z}_i = (z_{i1}, ..., z_{iK})^T\),满足 \(z_{ik} \in \{0, 1\}\) 且 \(\sum_{k=1}^K z_{ik} = 1\)。当 \(z_{ik} = 1\) 时,表示数据点 \(\mathbf{x}_i\) 来自于第 \(k\) 个高斯成分。
隐变量 \(\mathbf{z}_i\) 的先验分布由混合权重决定:
\[p(z_{ik} = 1) = \pi_k \]
给定 \(z_{ik} = 1\) 时,\(\mathbf{x}_i\) 的条件分布为:
\[p(\mathbf{x}_i | z_{ik} = 1) = \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \]
则完全数据(观测数据 + 隐变量)的对数似然函数为:
\[\ln p(\mathbf{X}, \mathbf{Z} | \Theta) = \sum_{i=1}^{N} \sum_{k=1}^{K} z_{ik} \left[ \ln \pi_k + \ln \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \right] \]
由于 \(z_{ik}\) 未知,我们无法直接最大化此函数。
步骤3:EM算法的E步(期望步)
E步的目标是计算完全数据对数似然关于隐变量后验分布的期望,这个期望被称为Q函数(Q-function)。
- 计算在当前参数估计 \(\Theta^{\text{old}}\) 下,隐变量的后验概率(也称为“响应度”或“责任”),记作 \(\gamma(z_{ik})\):
\[\gamma(z_{ik}) = p(z_{ik}=1 | \mathbf{x}_i, \Theta^{\text{old}}) = \frac{\pi_k^{\text{old}} \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_k^{\text{old}}, \boldsymbol{\Sigma}_k^{\text{old}})}{\sum_{j=1}^{K} \pi_j^{\text{old}} \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_j^{\text{old}}, \boldsymbol{\Sigma}_j^{\text{old}})} \]
这个公式是贝叶斯定理的直接应用。\(\gamma(z_{ik})\) 量化了当前数据点 \(\mathbf{x}_i\) 由第 \(k\) 个高斯成分生成的概率。
- 计算Q函数,它是完全数据对数似然关于后验分布 \(p(\mathbf{Z} | \mathbf{X}, \Theta^{\text{old}})\) 的期望:
\[Q(\Theta, \Theta^{\text{old}}) = E_{\mathbf{Z} | \mathbf{X}, \Theta^{\text{old}}} [\ln p(\mathbf{X}, \mathbf{Z} | \Theta)] = \sum_{i=1}^{N} \sum_{k=1}^{K} \gamma(z_{ik}) \left[ \ln \pi_k + \ln \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k) \right] \]
注意,在求期望时,未知的 \(z_{ik}\) 被其期望 \(E[z_{ik}] = \gamma(z_{ik})\) 所替代。
步骤4:EM算法的M步(最大化步)
M步的目标是找到一组新的参数 \(\Theta^{\text{new}}\) 以最大化Q函数 \(Q(\Theta, \Theta^{\text{old}})\)。由于Q函数中关于 \(\pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k\) 的部分是分离的,我们可以分别对它们求导并令导数为零来求解。
- 更新混合权重 \(\pi_k\)
需要考虑约束 \(\sum_k \pi_k = 1\),使用拉格朗日乘子法。构建拉格朗日函数:
\[ \mathcal{L}_{\pi} = \sum_{i=1}^{N} \sum_{k=1}^{K} \gamma(z_{ik}) \ln \pi_k + \lambda (\sum_{k=1}^{K} \pi_k - 1) \]
对 \(\pi_k\) 求导并令其为零,再结合约束,解得:
\[ \pi_k^{\text{new}} = \frac{1}{N} \sum_{i=1}^{N} \gamma(z_{ik}) = \frac{N_k}{N} \]
其中 \(N_k = \sum_{i=1}^{N} \gamma(z_{ik})\) 可视为“属于”第 \(k\) 个成分的“有效”数据点数量。新权重正比于这个有效点数。
- 更新均值向量 \(\boldsymbol{\mu}_k\)
对Q函数中与 \(\boldsymbol{\mu}_k\) 相关的部分(即 \(\sum_i \gamma(z_{ik}) \ln \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k)\))关于 \(\boldsymbol{\mu}_k\) 求导,令导数为零向量:
\[ \frac{\partial}{\partial \boldsymbol{\mu}_k} \left[ -\frac{1}{2} \sum_{i=1}^{N} \gamma(z_{ik}) (\mathbf{x}_i - \boldsymbol{\mu}_k)^T \boldsymbol{\Sigma}_k^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k) \right] = 0 \]
解得:
\[ \boldsymbol{\mu}_k^{\text{new}} = \frac{1}{N_k} \sum_{i=1}^{N} \gamma(z_{ik}) \mathbf{x}_i \]
新均值是所有样本的加权平均,权重是各样本属于该成分的后验概率。
- 更新协方差矩阵 \(\boldsymbol{\Sigma}_k\)
同样对Q函数中与 \(\boldsymbol{\Sigma}_k\) 相关的部分求导(或利用多元高斯分布的MLE结果),得到:
\[ \boldsymbol{\Sigma}_k^{\text{new}} = \frac{1}{N_k} \sum_{i=1}^{N} \gamma(z_{ik}) (\mathbf{x}_i - \boldsymbol{\mu}_k^{\text{new}})(\mathbf{x}_i - \boldsymbol{\mu}_k^{\text{new}})^T \]
新协方差矩阵是样本与当前新均值之间偏差的加权外积的平均,权重同样是后验概率。
步骤5:迭代与收敛
将E步和M步反复迭代:
- 初始化模型参数 \(\Theta^{(0)} = \{ \pi_k^{(0)}, \boldsymbol{\mu}_k^{(0)}, \boldsymbol{\Sigma}_k^{(0)} \}\)。
- 重复以下步骤直到参数变化小于阈值或对数似然收敛:
a. E步:使用当前参数 \(\Theta^{(t)}\) 计算所有 \(\gamma(z_{ik})^{(t)}\)。
b. M步:利用计算出的 \(\gamma(z_{ik})^{(t)}\) 更新参数,得到 \(\Theta^{(t+1)}\)。 - 最终,我们得到模型参数的极大似然估计(或局部极大值点)。同时,隐变量的后验概率 \(\gamma(z_{ik})\) 给出了每个数据点属于各个成分的“软”分配,可以用于软聚类。
算法总结:EM算法通过迭代地“填补”缺失的隐变量信息(E步,计算期望)和“优化”模型参数(M步,最大化期望)来求解GMM的参数。它保证了每一步迭代都不会降低观测数据的对数似然值,并最终收敛到一个局部最优解。