期望最大化（EM）算法在混合t分布（Mixture of t-Distributions）中的参数估计过程
题目描述
混合t分布模型是一种用于对具有重尾特性（即存在较多异常值）的数据进行建模的概率模型。与高斯混合模型（GMM）假设每个成分服从高斯分布（轻尾）不同，混合t分布假设每个混合成分服从t分布。t分布的自由度参数控制着分布的尾部厚度：自由度越小，尾部越厚，对异常值越稳健。由于t分布可以看作是无限个方差不相同的高斯分布的混合，其参数估计无法直接通过解析方法完成。
期望最大化（EM）算法是解决此类含隐变量模型参数估计的经典方法。本题将详细讲解如何使用EM算法估计混合t分布中的参数，包括：混合权重、每个成分的均值向量、协方差矩阵以及自由度参数。
解题过程（逐步讲解）
1. 模型设定与符号定义
假设我们有观测数据集：\(\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N\}\)，其中每个\(\mathbf{x}_i \in \mathbb{R}^D\)。
- 混合成分数：\(K\)
- 隐变量：
- 成分指示变量\(z_i \in \{1,\dots,K\}\)，表示第\(i\)个样本来自哪个成分。\(P(z_i = k) = \pi_k\)，满足\(\sum_{k=1}^K \pi_k = 1, \pi_k \ge 0\)。
- 辅助变量\(u_{ik}\)（来自t分布的表示）：t分布可以表示为一个高斯分布与其精度（方差的倒数）的混合，其中精度服从Gamma分布。即：
\[ \mathbf{x}_i | z_i=k, u_{ik} \sim \mathcal{N}(\boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k / u_{ik}), \quad u_{ik} | z_i=k \sim \text{Gamma}(\nu_k/2, \nu_k/2) \]
其中$\nu_k > 0$是t分布的自由度，控制尾部厚度。
- 待估参数：\(\boldsymbol{\theta} = \{\pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k, \nu_k\}_{k=1}^K\)。
2. 完全数据似然函数
给定隐变量\(z_i\)和\(u_{ik}\)，观测\(\mathbf{x}_i\)的条件分布是高斯分布。完全数据的联合分布为：
\[p(\mathbf{x}_i, u_{ik}, z_i=k | \boldsymbol{\theta}) = \pi_k \cdot \mathcal{N}(\mathbf{x}_i | \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k/u_{ik}) \cdot \text{Gamma}(u_{ik} | \nu_k/2, \nu_k/2) \]
其中Gamma分布的密度函数为：\(\text{Gamma}(u | a, b) = \frac{b^a}{\Gamma(a)} u^{a-1} e^{-b u}\)。
3. EM算法框架回顾
EM算法通过迭代以下两步求解：
- E步：在给定当前参数\(\boldsymbol{\theta}^{\text{old}}\)下，计算完全数据的对数似然关于隐变量的条件期望（即Q函数）：
\[ Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text{old}}) = \mathbb{E}_{p(\mathbf{z},\mathbf{u}|\mathbf{X},\boldsymbol{\theta}^{\text{old}})}[\log p(\mathbf{X},\mathbf{z},\mathbf{u}|\boldsymbol{\theta})] \]
- M步：最大化Q函数，更新参数：
\[ \boldsymbol{\theta}^{\text{new}} = \arg\max_{\boldsymbol{\theta}} Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text{old}}) \]
4. E步：计算需要的期望量
我们需要计算以下两个期望：
- 成分指示变量的后验概率（责任）：
\[ \gamma_{ik} = P(z_i = k | \mathbf{x}_i, \boldsymbol{\theta}^{\text{old}}) \]
由贝叶斯公式：
\[ \gamma_{ik} = \frac{\pi_k^{\text{old}} \cdot t(\mathbf{x}_i | \boldsymbol{\mu}_k^{\text{old}}, \boldsymbol{\Sigma}_k^{\text{old}}, \nu_k^{\text{old}})}{\sum_{j=1}^K \pi_j^{\text{old}} \cdot t(\mathbf{x}_i | \boldsymbol{\mu}_j^{\text{old}}, \boldsymbol{\Sigma}_j^{\text{old}}, \nu_j^{\text{old}})} \]
其中\(t(\cdot)\)是多元t分布密度函数：
\[ t(\mathbf{x} | \boldsymbol{\mu}, \boldsymbol{\Sigma}, \nu) = \frac{\Gamma\left(\frac{\nu + D}{2}\right)}{(\pi \nu)^{D/2} \Gamma(\nu/2) |\boldsymbol{\Sigma}|^{1/2}} \left[1 + \frac{1}{\nu} (\mathbf{x} - \boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) \right]^{-(\nu+D)/2} \]
- 辅助变量\(u_{ik}\)的条件期望：
给定\(z_i=k\)和\(\mathbf{x}_i\)，\(u_{ik}\)的后验分布仍为Gamma分布：
\[ p(u_{ik} | \mathbf{x}_i, z_i=k, \boldsymbol{\theta}^{\text{old}}) = \text{Gamma}\left(u_{ik} \left| \frac{\nu_k^{\text{old}} + D}{2}, \frac{\nu_k^{\text{old}} + \delta_{ik}^{\text{old}}}{2} \right. \right) \]
其中\(\delta_{ik}^{\text{old}} = (\mathbf{x}_i - \boldsymbol{\mu}_k^{\text{old}})^T (\boldsymbol{\Sigma}_k^{\text{old}})^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k^{\text{old}})\)是马氏距离的平方。
Gamma分布的期望为形状参数除以率参数，因此：
\[ \hat{u}_{ik} = \mathbb{E}[u_{ik} | \mathbf{x}_i, z_i=k, \boldsymbol{\theta}^{\text{old}}] = \frac{\nu_k^{\text{old}} + D}{\nu_k^{\text{old}} + \delta_{ik}^{\text{old}}} \]
另外，我们还需要\(\mathbb{E}[\log u_{ik} | \cdot]\)，它涉及Digamma函数\(\psi(\cdot)\)：
\[ \hat{\eta}_{ik} = \mathbb{E}[\log u_{ik} | \mathbf{x}_i, z_i=k, \boldsymbol{\theta}^{\text{old}}] = \psi\left( \frac{\nu_k^{\text{old}} + D}{2} \right) - \log\left( \frac{\nu_k^{\text{old}} + \delta_{ik}^{\text{old}}}{2} \right) \]
5. M步：更新参数
最大化Q函数得到各参数的更新公式：
- 混合权重\(\pi_k\)：
\[ \pi_k^{\text{new}} = \frac{1}{N} \sum_{i=1}^N \gamma_{ik} \]
- 均值向量\(\boldsymbol{\mu}_k\)：
\[ \boldsymbol{\mu}_k^{\text{new}} = \frac{\sum_{i=1}^N \gamma_{ik} \hat{u}_{ik} \mathbf{x}_i}{\sum_{i=1}^N \gamma_{ik} \hat{u}_{ik}} \]
注意：这里每个样本以权重\(\gamma_{ik} \hat{u}_{ik}\)参与加权平均，异常值（马氏距离大）的\(\hat{u}_{ik}\)较小，因此对均值估计的影响被削弱，体现了t分布的稳健性。
- 协方差矩阵\(\boldsymbol{\Sigma}_k\)：
\[ \boldsymbol{\Sigma}_k^{\text{new}} = \frac{\sum_{i=1}^N \gamma_{ik} \hat{u}_{ik} (\mathbf{x}_i - \boldsymbol{\mu}_k^{\text{new}})(\mathbf{x}_i - \boldsymbol{\mu}_k^{\text{new}})^T}{\sum_{i=1}^N \gamma_{ik}} \]
同样，协方差估计也被权重\(\hat{u}_{ik}\)稳健化。
- 自由度\(\nu_k\)：
\(\nu_k\)的更新没有封闭解，需要通过求解以下方程得到：
\[ \log\left(\frac{\nu_k}{2}\right) + 1 - \psi\left(\frac{\nu_k}{2}\right) + \frac{1}{\sum_{i=1}^N \gamma_{ik}} \sum_{i=1}^N \gamma_{ik} \left( \hat{\eta}_{ik} - \hat{u}_{ik} \right) = 0 \]
其中\(\psi\)是Digamma函数。这个方程通常用数值方法（如牛顿法）求解。注意，\(\nu_k\)控制尾部厚度：若数据在该成分中异常值多，则估计的\(\nu_k\)会较小（厚尾）；若数据近似高斯，则\(\nu_k\)会较大。
6. 算法流程总结
- 初始化：随机或通过K-means初始化\(\{\pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k\}\)，设\(\nu_k\)为较大值（如20）。
- E步：
- 计算马氏距离\(\delta_{ik} = (\mathbf{x}_i - \boldsymbol{\mu}_k)^T \boldsymbol{\Sigma}_k^{-1} (\mathbf{x}_i - \boldsymbol{\mu}_k)\)。
- 计算责任\(\gamma_{ik}\)（用多元t密度公式）。
- 计算期望\(\hat{u}_{ik} = \frac{\nu_k + D}{\nu_k + \delta_{ik}}\)和\(\hat{\eta}_{ik} = \psi\left( \frac{\nu_k + D}{2} \right) - \log\left( \frac{\nu_k + \delta_{ik}}{2} \right)\)。
- M步：
- 更新\(\pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k\)使用上述封闭解。
- 对每个\(k\)，数值求解关于\(\nu_k\)的方程，更新自由度。
- 重复E步和M步，直到对数似然\(\log p(\mathbf{X}|\boldsymbol{\theta})\)的变化小于阈值。
核心要点
- t分布的稳健性：通过辅助变量\(u_{ik}\)实现，异常值的\(u_{ik}\)较小，在估计均值、协方差时权重降低。
- 自由度参数的意义：\(\nu_k\)是估计的关键，小\(\nu_k\)对应厚尾，能更好地容纳离群点。
- 与GMM的EM对比：GMM的EM是混合t分布EM在\(\nu_k \to \infty\)时的特例（此时\(\hat{u}_{ik} \to 1\)）。
- 计算复杂度：主要在于协方差矩阵求逆（计算\(\delta_{ik}\)）和求解\(\nu_k\)的方程，每次迭代为\(O(N K D^2)\)。
该算法广泛应用于金融数据、图像处理等领域，其中数据常包含异常值或呈现重尾特性。