前面介绍的聚类和话题分析方法都有具体的应用场景，而这一讲我们会转向更加基础的机器学习理论——含有隐变量的概率模型估计。

• 具体而言，概率模型估计方法主要用于解决无监督概率模型（如GMM、PLSA、LDA等）的学习问题，包括EM算法、变分推理和马尔可夫链蒙特卡罗法等。

EM算法

在之前的人工智能导论与数理统计Cheat Sheet中，我们都简单提到了EM（期望-最大化）算法。它作为迭代算法，用于含有隐变量的概率模型参数的极大似然估计或极大后验概率估计。下面我们将更深入地研究其原理（证明其收敛性），并给出其应用，以及变分EM算法。

• 首先，回顾一下EM算法的经典引例——三硬币投掷（背景可参见人工智能导论）。设模型参数$\boldsymbol{\theta}=(\pi,p,q)$（分别表示三个硬币$A,B,C$正面出现的概率），观测数据为$\mathbf{x}=x_1,\cdots,x_n$（$n$次投掷的结果），未观测数据为$\mathbf{z}=z_1,\cdots,z_n$（表示硬币$A$的投掷结果）。则观测数据的似然函数可表示为 $$\begin{aligned} P(\mathbf{x}|\boldsymbol{\theta})=\prod_{i=1}^nP(x_i|\boldsymbol{\theta})&=\prod_{i=1}^n\left(\sum_{z_i}P(z_i|\boldsymbol{\theta})P(x_i|z_i,\boldsymbol{\theta})\right)\\ &=\prod_{i=1}^n\left[\pi(p^{x_i}(1-p)^{(1-x_i)})+(1-\pi)q^{x_i}(1-q)^{x_i}\right] \end{aligned}$$ 则模型参数的极大似然估计为$\displaystyle\hat{\boldsymbol{\theta}}=\argmax_{\boldsymbol{\theta}}\log P(\mathbf{x}|\boldsymbol{\theta})$。其使用EM迭代算法求解，步骤如下：
  1. 选取参数初值$\boldsymbol{\theta}^{(0)}=(\pi^{(0)},p^{(0)},q^{(0)})$；
  2. 记第$t$次迭代参数估计值为$\boldsymbol{\theta}^{(t)}=(\pi^{(t)},p^{(t)},q^{(t)})$，那么在第$t+1$次迭代中：
     • E步：计算当前参数估计值下观测数据$x_i$来自硬币$B$的概率 $$\gamma_i^{(t)} = \frac{\pi^{(t)} (p^{(t)})^{x_i} (1 - p^{(t)})^{1 - x_i}}{\pi^{(t)} (p^{(t)})^{x_i} (1 - p^{(t)})^{1 - x_i} + (1 - \pi^{(t)}) (q^{(t)})^{x_i} (1 - q^{(t)})^{1 - x_i}}, \quad i = 1, 2, \ldots, n$$ 而来自硬币$C$的概率则为$1-\gamma_i^{(t)}$，进而得到未观测数据$z_i$的概率估计值。
     • M步：计算参数新估计值 $$\begin{aligned} \pi^{(t+1)} &= \frac{1}{n} \sum_{i=1}^{n} \gamma_i^{(t)} \\ p^{(t+1)} &= \frac{\sum_{i=1}^{n} \gamma_i^{(t)} x_i}{\sum_{i=1}^{n} \gamma_i^{(t)}} \\ q^{(t+1)} &= \frac{\sum_{i=1}^{n} (1 - \gamma_i^{(t)}) x_i}{\sum_{i=1}^{n} (1 - \gamma_i^{(t)})} \end{aligned}$$
  3. 由此不断迭代，直到参数收敛。

  注：EM算法与初值的选择有关，选择不同的初值可能得到不同的参数估计值。

• 下面再给出EM算法的标准形式（与数理统计Cheat Sheet中符号略有不同）：
  • 输入: 观测数据$\mathbf{x}$。
  • 输出: 模型参数$\hat{\boldsymbol{\theta}}$。
  1. 选择初值: 选择参数的初值$\boldsymbol{\theta}^{(0)}$，开始迭代。
  2. E 步: 记$\boldsymbol{\theta}^{(t)}$为第$t$次迭代参数$\theta$的估计值，在第$t+1$次迭代的E步，计算 $$\begin{aligned} Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)}) &= \mathbb{E}_{P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)})} [\log P(\mathbf{x}, \mathbf{z}|\boldsymbol{\theta})] \\ &= \sum_{z} P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) \log P(\mathbf{x}, \mathbf{z}|\boldsymbol{\theta}) \end{aligned}$$ 这里，$P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)})$是给定观测数据$\mathbf{x}$和参数$\boldsymbol{\theta}^{(t)}$条件下未观测数据$\mathbf{z}$的条件概率分布。
  3. M 步: 求使$Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})$最大化的$\theta$，确定第$t+1$次迭代的参数的估计值$\boldsymbol{\theta}^{(t+1)}$。 $$\boldsymbol{\theta}^{(t+1)} = \argmax_{\boldsymbol{\theta}}Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})$$
  4. 迭代判断: 重复步骤 2 和步骤 3，直到收敛。收敛条件通常是对较小的正数$\epsilon_1, \epsilon_2$， $$\|\boldsymbol{\theta}^{(t+1)} - \boldsymbol{\theta}^{(t)}\| < \epsilon_1 \quad \text{或} \quad \|Q(\boldsymbol{\theta}^{(t+1)}, \boldsymbol{\theta}^{(t)}) - Q(\boldsymbol{\theta}^{(t)}, \boldsymbol{\theta}^{(t)})\| < \epsilon_2$$
  5. 最终得到模型参数$\hat{\boldsymbol{\theta}} = \boldsymbol{\theta}^{(t+1)}$。
• 由此可见，EM算法的核心就是计算$Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})$，即$Q$函数，其本质是完全数据的对数似然函数$\log P(\mathbf{x}, \mathbf{z}|\boldsymbol{\theta})$关于未观测数据的后验概率分布$P(\mathbf{z}|\mathbf{x},\boldsymbol{\theta}^{(t)})$的期望。

算法原理

那么上述$Q$函数的形式是如何得到的？下面我们从似然函数开始推导：

• 考虑不完全数据$\mathbf{x}$关于$\boldsymbol{\theta}$的对数似然函数 $$L(\boldsymbol{\theta})=\log P(\mathbf{x}|\boldsymbol{\theta})=\log\sum_{\mathbf{z}}P(\mathbf{x}, \mathbf{z}|\boldsymbol{\theta})=\log\left(\sum_{\mathbf{z}}P(\mathbf{x}|\mathbf{z},\boldsymbol{\theta})P(\mathbf{z}|\boldsymbol{\theta})\right)$$ 其无法直接通过解析方法优化。对此，EM算法通过迭代更新$L(\boldsymbol{\theta})$的下界实现。
• 具体而言，设第$t$次迭代后似然函数为$L(\boldsymbol{\theta}^{(t)})$，我们希望更新参数后似然函数$L(\boldsymbol{\theta})\geq L(\boldsymbol{\theta}^{(t)})$，于是考虑二者之差： $$\begin{aligned} L(\boldsymbol{\theta}) - L(\boldsymbol{\theta}^{(t)}) &=\log\left(\sum_{\mathbf{z}}P(\mathbf{x}|\mathbf{z}, \boldsymbol{\theta}) P(\mathbf{z}|\boldsymbol{\theta})\right) - \log P(\mathbf{x}|\boldsymbol{\theta}^{(t)})\\ &= \log\left(\sum_{\mathbf{z}} P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) \frac{P(\mathbf{x}|\mathbf{z}, \boldsymbol{\theta}) P(\mathbf{z}|\boldsymbol{\theta})}{P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)})}\right) - \log P(\mathbf{x}|\boldsymbol{\theta}^{(t)}) \\ &\geq \sum_{\mathbf{z}} P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) \log \frac{P(\mathbf{x}|\mathbf{z}, \boldsymbol{\theta}) P(\mathbf{z}|\boldsymbol{\theta})}{P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)})} - \log P(\mathbf{x}|\boldsymbol{\theta}^{(t)}) \\ &= \sum_{\mathbf{z}} P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) \log \frac{P(\mathbf{x}|\mathbf{z}, \boldsymbol{\theta}) P(\mathbf{z}|\boldsymbol{\theta})}{P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) P(\mathbf{x}|\boldsymbol{\theta}^{(t)})} \end{aligned}$$ 其中不等号使用的是Jensen不等式。由此得到$L(\boldsymbol{\theta})$的一个下界 $$B(\boldsymbol{\theta},\boldsymbol{\theta}^{(t)})\triangleq L(\boldsymbol{\theta}^{(t)})+\sum_{\mathbf{z}} P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) \log \frac{P(\mathbf{x}|\mathbf{z}, \boldsymbol{\theta}) P(\mathbf{z}|\boldsymbol{\theta})}{P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) P(\mathbf{x}|\boldsymbol{\theta}^{(t)})}$$ 由此可知$B(\boldsymbol{\theta}^{(t)},\boldsymbol{\theta}^{(t)})=L(\boldsymbol{\theta}^{(t)})$，且我们可以通过极大化下界得到迭代参数，即 $$\begin{aligned} \boldsymbol{\theta}^{(t+1)}&=\argmax_{\boldsymbol{\theta}}B(\boldsymbol{\theta},\boldsymbol{\theta}^{(t)})\\ &= \argmax_{\boldsymbol{\theta}} \left( L(\boldsymbol{\theta}^{(t)}) + \sum_{\mathbf{z}} P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) \log \frac{P(\mathbf{x}|\mathbf{z}, \boldsymbol{\theta}) P(\mathbf{z}|\boldsymbol{\theta})}{P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) P(\mathbf{x}|\boldsymbol{\theta}^{(t)})} \right) \\ &= \argmax_{\boldsymbol{\theta}} \left( \sum_{\mathbf{z}} P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) \log \left( P(\mathbf{x}|\mathbf{z}, \boldsymbol{\theta}) P(\mathbf{z}|\boldsymbol{\theta}) \right) \right) \\ &= \argmax_{\boldsymbol{\theta}} \left( \sum_{\mathbf{z}} P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) \log P(\mathbf{x}, \mathbf{z}|\boldsymbol{\theta}) \right) \\ &= \argmax_{\boldsymbol{\theta}} Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)}) \end{aligned}$$ 这样就得到了EM算法的M步，同时也解释了为什么算法对初值敏感（局部最优不一定全局最优）。

算法收敛性

• EM算法的收敛性包括对数似然函数$L(\boldsymbol{\theta}^{(t)})$的收敛性和参数$\boldsymbol{\theta}^{(t)}$本身的收敛性。下面我们主要讨论前者的收敛性证明（关于参数收敛性需要增加一定的条件，可参见吴建福的论文）
• 由于$L(\boldsymbol{\theta}^{(t)})=\log P(\mathbf{x}|\boldsymbol{\theta}^{(t)})$显然有上界，因此根据单调有界定理，只需要证明其单调递增，即 $$\log P(\mathbf{x}|\boldsymbol{\theta}^{(t+1)})\geq \log P(\mathbf{x}|\boldsymbol{\theta}^{(t)})$$ 其证明如下：
  • 因为$P(\mathbf{x}|\boldsymbol{\theta})=\dfrac{P(\mathbf{x},\mathbf{z}|\boldsymbol{\theta})}{P(\mathbf{z}|\mathbf{x},\boldsymbol{\theta})}$，且 $$Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})= \sum_{z} P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) \log P(\mathbf{x}, \mathbf{z}|\boldsymbol{\theta})$$ 定义 $$H(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})\triangleq \sum_{z} P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) \log P(\mathbf{z}|\mathbf{x},\boldsymbol{\theta})$$ 那么有$\log P(\mathbf{x}|\boldsymbol{\theta})=Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})-H(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})$，于是 $$\log P(\mathbf{x}|\boldsymbol{\theta}^{(t+1)}) - \log P(\mathbf{x}|\boldsymbol{\theta}^{(t)})= \left[ Q(\boldsymbol{\theta}^{(t+1)}, \boldsymbol{\theta}^{(t)}) - Q(\boldsymbol{\theta}^{(t)}, \boldsymbol{\theta}^{(t)}) \right] - \left[ H(\boldsymbol{\theta}^{(t+1)}, \boldsymbol{\theta}^{(t)}) - H(\boldsymbol{\theta}^{(t)}, \boldsymbol{\theta}^{(t)}) \right]$$
  • 由之前的原理推导可知$Q$函数在迭代过程中单调递增，而对于$H$函数： $$\begin{aligned} H(\boldsymbol{\theta}^{(t+1)}, \boldsymbol{\theta}^{(t)}) - H(\boldsymbol{\theta}^{(t)}, \boldsymbol{\theta}^{(t)}) &= \sum_{\mathbf{z}} P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) \log \frac{P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t+1)})}{P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)})} \\ &\leq \log \left( \sum_{\mathbf{z}} P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) \cdot \frac{P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t+1)})}{P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)})} \right) \\ &= \log \left( \sum_{\mathbf{z}} P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t+1)}) \right) = \log 1 = 0 \end{aligned}$$ 即$H$函数在迭代过程中单调递减。综上可得似然函数单调递增，从而收敛。

广义EM算法

• 注意到，上述EM算法在处理隐变量时直接使用后验概率分布$P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)})$，但现实情况中这个后验分布往往不容易直接求出。【和上述算法原理中似然函数无法直接解析优化原因类似】于是，我们将EM算法进行推广，让后验分布也加入迭代过程中，这便是广义EM算法（GEM）。
• 在介绍广义EM算法之前，先引入$F$函数：假设未观测数据$\mathbf{z}$的概率分布为$P(\mathbf{z}|\boldsymbol{\theta})$，定义分布$\tilde{P}$与参数$\boldsymbol{\theta}$的函数$F(\tilde{P}, \boldsymbol{\theta})$如下： $$F(\tilde{P}, \boldsymbol{\theta}) = \mathbb{E}_{\tilde{P}(\mathbf{z}|\boldsymbol{\theta})}[\log P(\mathbf{x}, \mathbf{z}|\boldsymbol{\theta})] + H(\tilde{P})$$ 其中$H(\tilde{P}) = -\mathbb{E}_{\tilde{P}(\mathbf{z}|\boldsymbol{\theta})} \log \tilde{P}(\mathbf{z}|\boldsymbol{\theta})$是分布$\tilde{P}(\mathbf{z}|\boldsymbol{\theta})$的信息熵。
• $F$函数具有以下性质：对于固定的$\boldsymbol{\theta}$，若存在分布$\tilde{P}$使$F(\tilde{P},\boldsymbol{\theta})$最大，则一定有 $$\tilde{P}(\mathbf{z}|\boldsymbol{\theta}) = P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta})$$ 此时$F(\tilde{P}, \boldsymbol{\theta})=\log P(\mathbf{x}|\boldsymbol{\theta})$。
  • 实际上，$F$函数还可以等价写为对数似然与KL散度之和，推导如下： $$\begin{aligned} F(\tilde{P}, \boldsymbol{\theta}) &= \mathbb{E}_{\tilde{P}(\mathbf{z}|\boldsymbol{\theta})}[\log P(\mathbf{x}, \mathbf{z}|\boldsymbol{\theta})] + H(\tilde{P}) \\ &= \mathbb{E}_{\tilde{P}}[\log P(\mathbf{x}, \mathbf{z}|\boldsymbol{\theta})] - \mathbb{E}_{\tilde{P}}[\log \tilde{P}(\mathbf{z}|\boldsymbol{\theta})] \\ &= \mathbb{E}_{\tilde{P}}[\log P(\mathbf{x}|\boldsymbol{\theta})] + \mathbb{E}_{\tilde{P}}[\log P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta})] - \mathbb{E}_{\tilde{P}}[\log \tilde{P}(\mathbf{z}|\boldsymbol{\theta})] \\ &= \log P(\mathbf{x}|\boldsymbol{\theta}) - \left( \mathbb{E}_{\tilde{P}}[\log \tilde{P}(\mathbf{z}|\boldsymbol{\theta})] - \mathbb{E}_{\tilde{P}}[\log P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta})] \right) \\ &= \log P(\mathbf{x}|\boldsymbol{\theta}) - \mathbb{E}_{\tilde{P}}\left[ \log \frac{\tilde{P}(\mathbf{z}|\boldsymbol{\theta})}{P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta})} \right] \\ &= \log P(\mathbf{x}|\boldsymbol{\theta}) - \text{KL}\left(\tilde{P}(\mathbf{z}|\boldsymbol{\theta}) \parallel P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta})\right) \end{aligned}$$ 这也能解释上述性质为何成立。【KL散度非负，当且仅当两个分布相同时取$0$】

最大化-最大化算法

有了$F$函数的定义，我们就能将EM算法扩展为最大化-最大化算法。具体形式如下（本质是交替迭代）：

• 初始化参数$\boldsymbol{\theta}^{(0)}$；

• 设$\boldsymbol{\theta}^{(t)}$为第$t$次迭代参数$\boldsymbol{\theta}$的估计，$\tilde{P}^{(t)}$为第$t$次迭代概率分布$\tilde{P}$的估计。第$t+1$次迭代的两步如下：

  1. 对固定的$\boldsymbol{\theta}^{(t)}$，求$\tilde{P}^{(t+1)}$使$F(\tilde{P}, \boldsymbol{\theta}^{(t)})$最大化。
     由上述$F$函数性质可知，当$F(\tilde{P}, \boldsymbol{\theta}^{(t)})$最大时 $$\tilde{P}^{(t+1)}\triangleq\tilde{P}(\mathbf{z}|\boldsymbol{\theta}^{(t)}) = P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)})$$ 因此 $$\begin{aligned} F(\tilde{P}^{(t+1)}, \boldsymbol{\theta}) &= \mathbb{E}_{\tilde{P}^{(t+1)}}[\log P(\mathbf{x}, \mathbf{z}|\boldsymbol{\theta})] + H(\tilde{P}^{(t+1)}) \\ &= \sum_{\mathbf{z}} P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)}) \log P(\mathbf{x}, \mathbf{z}|\boldsymbol{\theta}) + H(\tilde{P}^{(t+1)})\\ &= Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(t)})+H(\tilde{P}^{(t+1)}) \end{aligned}$$
  2. 对固定的$\tilde{P}^{(t+1)}$，求$\boldsymbol{\theta}^{(t+1)}$使$F(\tilde{P}^{(t+1)}, \boldsymbol{\theta})$最大化。
     由上述推导可知 $$\boldsymbol{\theta}^{(t+1)}=\argmax_{\boldsymbol{\theta}}F(\tilde{P}^{(t+1)},\boldsymbol{\theta})=\argmax_{\boldsymbol{\theta}}Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(t)})$$

  由此可知，由EM算法与$F$函数的最大化-最大化算法得到的参数估计序列是一致的。

• 重复上述迭代过程，直到收敛，就得到参数估计值$\hat{\boldsymbol{\theta}}$。

EM算法应用：高斯混合模型

• 下面给出一个EM算法的经典应用：高斯混合模型（Gaussian Mixture Model，GMM）的学习。模型的本质是多个高斯分布（即正态分布）模型的线性组合，表达式如下：

  $$p(x|\boldsymbol{\theta})=\sum_{j=1}^k\pi_jf(x|\boldsymbol{\theta}_j)$$

  其中$\pi_1,\cdots,\pi_k$构成一个概率分布，$f(x|\boldsymbol{\theta}_j)$为第$j$个正态分布（模型）的密度函数，参数为$\boldsymbol{\theta}_j=(\mu_j,\sigma^2_j)$。【$k$属于超参数】

• 高斯混合模型属于生成模型，主要用于数据聚类，与k均值聚类为代表的硬聚类不同，其属于“软聚类”，没有严格的划分边界。

• 接下来我们将用EM算法估计模型的参数$\boldsymbol{\theta}$，包括$(\pi_j,\mu_j,\sigma_j^2)_{j=1}^k$。具体步骤如下：

  1. 求完全数据的对数似然函数
     设$x_1,\cdots,x_n$为观测数据，引入隐变量$z_{ij}$，表示第$i$个观测值是否属于第$j$个高斯模型（如果属于取值$1$，否则取值$0$）。【满足$\displaystyle\sum_{j=1}^kz_{ij}=1$】
     • 考虑完全数据$(x_i,z_{i1},\cdots,z_{ij}),i=1,2,\cdots,n$，可写出其似然函数 $$\begin{aligned} p(\mathbf{x}, \mathbf{z}|\boldsymbol{\theta}) &= \prod_{i=1}^n p(x_i, z_{i1}, z_{i2}, \cdots, z_{ij}|\boldsymbol{\theta}) \\ &= \prod_{i=1}^n \prod_{j=1}^k \bigl( \pi_j f(x_i|\boldsymbol{\theta}_j) \bigr)^{z_{ij}} \\ &= \prod_{j=1}^k \prod_{i=1}^n \pi_j^{z_{ij}}\bigl( f(x_i|\boldsymbol{\theta}_j) \bigr)^{z_{ij}} \\ &= \prod_{j=1}^k \prod_{i=1}^n \pi_j^{z_{ij}} \left\{ \frac{1}{\sqrt{2\pi}\sigma_j} \exp\left[ -\frac{(x_i - \mu_j)^2}{2\sigma_j^2} \right] \right\}^{z_{ij}} \end{aligned}$$ 于是对数似然函数为 $$\log p(\mathbf{x}, \mathbf{z}|\boldsymbol{\theta}) = \sum_{j=1}^k \left\{ \sum_{i=1}^n z_{ij} \log \pi_j + \sum_{i=1}^n z_{ij} \left[ \log \frac{1}{\sqrt{2\pi}} - \log \sigma_j - \frac{1}{2\sigma_j^2}(\mathbf{x}_i - \mu_j)^2 \right] \right\}$$
  2. E步：求$Q$函数
     将上述对数似然函数对隐变量分布求期望得 $$\begin{aligned} Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)}) &= \mathbb{E}_{p(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)})} \left( \sum_{j=1}^{k} \left\{ \sum_{i=1}^{n} z_{ij} \log \pi_j + \sum_{i=1}^{n} z_{ij} \left[ \log \frac{1}{\sqrt{2\pi}} - \log \sigma_j - \frac{1}{2\sigma_j^2} (\mathbf{x}_i - \mu_j)^2 \right] \right\} \right) \\ &= \sum_{j=1}^{k} \left\{ \sum_{i=1}^{n} \mathbb{E}_{p(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)})}(z_{ij}) \log \pi_j + \sum_{i=1}^{n} \mathbb{E}_{p(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)})}(z_{ij}) \left[ \log \frac{1}{\sqrt{2\pi}} - \log \sigma_j - \frac{1}{2\sigma_j^2} (\mathbf{x}_i - \mu_j)^2 \right] \right\} \end{aligned}$$ 其中 $$\begin{aligned} \mathbb{E}_{p(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)})}(z_{ij}) &= p(z_{ij} = 1|\mathbf{x}, \boldsymbol{\theta}^{(t)}) \\ &= \frac{p(z_{ij} = 1, \mathbf{x}_i|\boldsymbol{\theta}^{(t)})} {\displaystyle\sum_{j=1}^k p(z_{ij} = 1, \mathbf{x}_i|\boldsymbol{\theta}^{(t)})} \\ &= \frac{p(z_{ij} = 1|\boldsymbol{\theta}^{(t)})p(\mathbf{x}_i|z_{ij} = 1, \boldsymbol{\theta}^{(t)})} {\displaystyle\sum_{j=1}^k p(z_{ij} = 1|\boldsymbol{\theta}^{(t)})p(\mathbf{x}_i|z_{ij} = 1, \boldsymbol{\theta}^{(t)})} \\ &= \frac{\pi_j^{(t)}f(\mathbf{x}_i|\boldsymbol{\theta}_j^{(t)})} {\displaystyle\sum_{j=1}^k \pi_j^{(t)}f(\mathbf{x}_i|\boldsymbol{\theta}_j^{(t)})}, \qquad i = 1, 2, \dots, n, \quad j = 1, 2, \dots, k \end{aligned}$$ 将其记作$\gamma_{ij}^{(t)}$，表示当前模型参数下第$i$个观测数据来自第$j$个分模型的后验概率，也被称作责任度（responsibility）。
  3. M步：$Q$函数最大化
     使用通式 $$\boldsymbol{\theta}^{(t+1)} = \argmax_{\boldsymbol{\theta}}Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})$$ 由于上述$Q$函数为凸函数，因此直接将其对各参数求偏导并使其为$0$即得参数新迭代值。具体结果如下： $$\begin{aligned} \mu_j^{(t+1)} &= \frac{\displaystyle\sum_{i=1}^n \gamma_{ij}^{(t)} \mathbf{x}_i}{\displaystyle\sum_{i=1}^n \gamma_{ij}^{(t)}}, \quad j = 1, 2, \dots, k \\[8pt] (\sigma_j^2)^{(t+1)} &= \frac{\displaystyle\sum_{i=1}^n \gamma_{ij}^{(t)} (\mathbf{x}_i - \mu_j^{(t+1)})^2}{\displaystyle\sum_{i=1}^n \gamma_{ij}^{(t)}}, \quad j = 1, 2, \dots, k \\[8pt] \pi_j^{(t+1)} &= \frac{\displaystyle\sum_{i=1}^n \gamma_{ij}^{(t)}}{n}, \quad j = 1, 2, \dots, k \end{aligned}$$

  重复上述迭代直到参数收敛即可。

• 实际上，当上述高斯混合过程中的EM算法里所有方差参数$\sigma_j^2=\epsilon\to 0$，那么这一算法就退化为了k均值聚类算法。

变分推理

前面提到的广义EM算法中，虽然没有直接求隐变量的后验分布，但在迭代过程中仍然需要计算当前参数下的后验分布$P(\mathbf{z}|\mathbf{x}, \boldsymbol{\theta}^{(t)})$，这往往是难以直接实现的。

• 目前主流的两种解决策略如下：一种是限制后验分布的搜索空间（代表是变分推理），另一种是使用随机抽样近似计算（代表是MCMC）。这里我们先介绍前一种方法。

变分贝叶斯方法

• 变分贝叶斯方法的核心思想是：设定一个后验分布空间$\mathcal{Q}$，使用空间中的元素$q(\mathbf{z})$（也称为变分分布）来近似后验分布$p(\mathbf{z}|\mathbf{x})$。
• 二者的相似度使用KL散度$\text{KL}(q(\mathbf{z})\parallel p(\mathbf{z}|\mathbf{x}))$计算，具体如下： $$\begin{aligned} \text{KL}(q(\mathbf{z}) \parallel p(\mathbf{z}|\mathbf{x})) &= \mathbb{E}_{q(\mathbf{z})} \left[ \log q(\mathbf{z}) \right] - \mathbb{E}_{q(\mathbf{z})} \left[ \log p(\mathbf{z}|\mathbf{x}) \right] \\ &= \mathbb{E}_{q(\mathbf{z})} \left[ \log q(\mathbf{z}) \right] - \mathbb{E}_{q(\mathbf{z})} \left[ \log p(\mathbf{x}, \mathbf{z}) \right] + \log p(\mathbf{x}) \\ &= \log p(\mathbf{x}) - \left\{ \mathbb{E}_{q(\mathbf{z})} \left[ \log p(\mathbf{x}, \mathbf{z}) \right] - \mathbb{E}_{q(\mathbf{z})} \left[ \log q(\mathbf{z}) \right] \right\} \end{aligned}$$ 又因为KL散度非负，因此有 $$\log p(\mathbf{x})\geq \left\{ \mathbb{E}_{q(\mathbf{z})} \left[ \log p(\mathbf{x}, \mathbf{z}) \right] - \mathbb{E}_{q(\mathbf{z})} \left[ \log q(\mathbf{z}) \right] \right\}$$ 将上式左端称作证据（evidence），右端记作$L(q)$，称作证据下界（evidence lower bound, ELBO）。【实际上，ELBO的形式和上述GEM的$F$函数完全一致】
• 因为KL散度越小，近似效果越好，所以变分贝叶斯方法的目标是KL散度最小化，基于上述推导我们可以将其转化为ELBO的最大化。【证据为常量】
• 变分分布空间$\mathcal{Q}$一般取为平均场空间，即其中每个元素$q(\mathbf{z})$都满足（条件）独立假设，即 $$q(\mathbf{z})=\prod_{i=1}^nq(z_i)$$

变分EM算法

在上述变分贝叶斯方法中，一般使用迭代方法最大化ELBO，此时也称为变分EM算法，可看作EM算法的推广。

• 变分EM算法的具体流程如下：
  • 输入：观测数据。
  • 输出：模型参数。
  1. 初始化参数$\boldsymbol{\theta}^{(0)}$。
  2. E 步：固定参数$\boldsymbol{\theta}^{(t)}$，求使$L(q(\mathbf{z}), \boldsymbol{\theta}^{(t)})$最大化的变分分布$q(\mathbf{z})$，得到$q^{(t+1)}(\mathbf{z})$。
  3. M 步：固定变分分布$q^{(t+1)}(\mathbf{z})$，求使$L(q^{(t+1)}(\mathbf{z}), \boldsymbol{\theta})$最大化的参数$\boldsymbol{\theta}$，得到$\boldsymbol{\theta}^{(t+1)}$。
  4. 令$t = t + 1$，重复步骤 2 和步骤 3，直到收敛。
  5. 给出模型参数$\hat{\boldsymbol{\theta}}$。
• 可以证明上述$L(q(\mathbf{z}), \boldsymbol{\theta}^{(t)})$序列一定单调递增且收敛。【但不保证收敛到全局最优解】
• 变分EM算法与传统EM算法之间的区别可以理解为贝叶斯学派与频率学派对参数估计的不同理解。EM算法往往用于模型简单的情况，一般估计的准确率更高；而变分EM算法用于模型复杂的情况，近似更多。

马尔可夫链蒙特卡罗法（MCMC）

上述EM算法都是通过迭代的方式进行概率模型参数估计，而下面我们将介绍另一种估计方法——统计模拟。

• 在数理统计Cheat Sheet中，我们提到蒙特卡罗分布通过随机抽样对参数进行近似估计；而在随机过程Cheat Sheet中，我们提到可以构造马尔科夫链使其平稳分布能够近似任意离散概率分布。【即Metropolis-Hastings算法，这也是最基本的马尔可夫链蒙特卡罗法】
• 下面我们会对MCMC进行进一步阐述：

蒙特卡罗法（Monte Carlo method）

蒙特卡罗法是统计学与机器学习中广泛使用的估计方法，其核心是随机抽样（random sampling）。

• 蒙特卡罗法除了直接抽样法之外，还包括接受-拒绝抽样法、重要性抽样法等，其适用于概率密度函数复杂、无法直接抽样的情形。

• 以单变量接受-拒绝抽样法为例，其算法步骤如下：

  1. 设随机变量$x$的概率密度函数为$p(x)$，另外找一个可以直接抽样的分布$q(x)$（也称为建议分布），要求存在$c>0$使$c\cdot q(x)\geq p(x)$在支撑集上总成立；
  2. 按照建议分布$q(x)$随机抽样得到样本$x_i$，再按照均匀分布在$(0,1)$范围内抽样得到$u$；
  3. 若$u\leq\dfrac{p(x_i)}{c\cdot q(x_i)}$，则将$x_i$作为抽样结果；否则回到上一步；
  4. 重复上述步骤直到得到$n$个样本。

  接受-拒绝法的优点是容易实现，缺点是效率可能不高。

• 蒙特卡罗法也可以用于数学期望估计：设随机变量$x$的概率密度函数为$p(x)$，目标是求函数$f(x)$关于$p(x)$的数学期望$\mathbb{E}_{p(x)}[f(x)]$。

  • 可以使用直接抽样法：按照概率分布$p(x)$独立抽取$n$个样本$x_1,\cdots,x_n$，分别计算其函数值$f(x_i)$，再求其平均。根据大数定律，

    $$\frac{1}{n}\sum_{i=1}^nf(x_i)\stackrel{p}{\longrightarrow}\mathbb{E}_{p(x)}[f(x)].$$

  • 当$p(x)$较复杂时，可采用重要性抽样法，其具体步骤如下（与接受-拒绝抽样法类似）：

    1. 选择概率密度函数为$q(x)$的概率分布作为建议分布，保证$q(x)$与目标分布$p(x)$有相同支撑集；
    2. 按照建议分布$q(x)$随机抽样$n$个样本$x_i, \quad i = 1, 2, \cdots, n$。
    3. 对于每个样本计算其重要性： $$w_i = \frac{p(x_i)}{q(x_i)}$$
    4. 计算所有样本的函数值的加权平均，得到函数期望的估计： $$\hat{f}_n = \frac{1}{n} \sum_{i=1}^{n} w_i f(x_i)$$

    重要性抽样法的优点同样是容易实现，缺点是函数期望估计的方差可能很大。

• 蒙特卡罗法还可以用于定积分的近似计算，称为蒙特卡罗积分。实际上，对于下述函数积分

  $$\int_{\mathcal{X}}h(x)\mathrm{d}x$$

  如果$h(x)$可分解为一个概率密度函数$p(x)$（支撑集为$\mathcal{X}$）与另一函数$f(x)$的乘积，那么有

  $$\int_{\mathcal{X}}h(x)\mathrm{d}x=\int_{\mathcal{X}}p(x)f(x)\mathrm{d}x=\mathbb{E}_{p(x)}[f(x)].$$

  这样我们就将函数的积分转换为另一个函数的数学期望的形式，然后就能套用上述蒙特卡罗估计方法。

• 上述积分计算对于机器学习（尤其是贝叶斯学习）非常有用：假设观测数据由随机变量$d \in \mathcal{D}$表示，模型由随机变量$m \in \mathcal{M}$表示，贝叶斯学习通过贝叶斯定理计算给定数据条件下模型的后验概率，并选择后验概率最大的模型。后验概率为

  $$p(m|d) = \frac{p(m)p(d|m)}{\int_{\mathcal{M}} p(m')p(d|m')\mathrm{d}m'}$$

  计算中需要归一化计算：

  $$p(d) = \int_{\mathcal{M}} p(d|m')p(m')\mathrm{d}m'$$

  如果有隐变量$z \in \mathcal{Z}$，还需要边缘化计算：

  $$p(m|d) = \int_{\mathcal{Z}} p(m, z|d)\mathrm{d}z$$

  如果有一个函数$f(m)$，可以计算该函数关于后验概率分布的数学期望：

  $$\mathbb{E}_{p(m|d)}[f(m)] = \int_{\mathcal{M}} f(m)p(m|d)\mathrm{d}m$$

  当观测数据和模型都很复杂的时候，以上积分的计算都比较困难。使用蒙特卡洛法（一般为MCMC）可以有效解决这种复杂计算问题。

马尔科夫链

上述蒙特卡罗方法得到的抽样样本都是相互独立的，但往往无法解决处理复杂的概率分布（如随机变量是多元、密度函数是非标准形式、随机变量各分量不独立等）。对此，我们引入马尔科夫链，通过其平稳分布逼近实际概率分布。

• 关于离散状态马尔科夫链的相关概念可参见随机过程Cheat Sheet，下面只简单说明连续状态马尔科夫链：
  • 定义连续状态马尔科夫链$\mathbf{X}=X_0X_1\cdots X_t\cdots$，其中$X_t$定义在连续状态空间$\mathcal{S}$，其状态转移概率分布由概率转移核（transition kernel）表示。
  • 对任意的$x\in\mathcal{S},A\subset\mathcal{S}$，转移核$P(x,A)=P(X_t\in A|X_{t-1}=x)$定义为 $$P(x,A)=\int_Ap(x,y)\mathrm{d}y$$ 其中$p(x,\cdot)$为概率密度函数，满足$\displaystyle P(x,\mathcal{S})=\int_{\mathcal{S}}p(x,y)\mathrm{d}y=1$。【有时也称$p(x,\cdot)$为转移核】
  • 连续状态马尔科夫链的平稳分布$\pi(x)$满足以下条件： $$\pi(y) = \int_{\mathcal{S}} p(x, y)\pi(x)\mathrm{d}x, \quad \forall y \in \mathcal{S} \quad\text{或}\quad\pi(A) = \int_{\mathcal{S}} P(x, A)\pi(x)\mathrm{d}x, \quad \forall A \subset \mathcal{S}$$

基本想法

马尔科夫链蒙特卡罗法的基本想法如下：在随机变量$x$的状态空间$\mathcal{S}$上定义一个满足遍历性定理的马尔可夫链$\mathbf{X} = X_0 X_1 \cdots X_t \cdots$，使其平稳分布就是抽样的目标分布$p(x)$。然后在这个马尔可夫链上进行随机游走，每个时刻得到一个样本。

• 根据遍历性定理，当时间趋于无穷时，样本状态的频率分布趋近平稳分布，样本的函数均值趋近函数的数学期望。
• 所以，当时间足够长（设$n,m$为正整数且足够大）时，随机游走得到的样本集合$\{x_{m+1}, x_{m+2}, \cdots, x_n\}$就是目标概率分布的抽样结果，得到的函数均值（遍历均值）就是要计算的数学期望值： $$\hat{f}_{m,n} = \frac{1}{n - m} \sum_{t=m+1}^{n} f(x_t)$$ 时刻$m$之前的时间段称为燃烧期（burn-in），$m$称为收敛步数（决定抽样无偏性），$n$称为迭代步数（决定均值计算精度）。
• MCMC收敛性的判断通常是经验性的，一般会每隔一段时间取一次样本，得到多个样本以后，计算遍历均值。当计算的均值稳定后，认为马尔可夫链已经收敛。也可以在马尔可夫链上并行进行多个随机游走，比较各个随机游走的遍历均值是否接近一致。
• MCMC得到样本序列具有相关性。在需要独立样本时，可以对其再次进行随机抽样，比如每隔一段时间取一次样本，将这样得到的子样本集合作为独立样本集合。
• MCMC比接受-拒绝法更容易实现，且其效率往往更高（样本拒绝率较低）。

常用的马尔可夫链蒙特卡罗法有Metropolis-Hastings算法与吉布斯抽样。下面进行详细介绍：

Metropolis-Hastings算法

Metropolis-Hastings算法已经在随机过程Cheat Sheet中简要提到过，下面我们对这一算法进行补充：

• 随机过程Cheat Sheet中的预选矩阵$\mathbf{Q}$其实就是建议分布，而$\alpha$矩阵则称为接受分布。在实际的马尔科夫链随机游走中，如果在时刻$t-1$处于状态$x$（即$\mathbf{X}_{t-1}=x$），则先按建议分布抽样产生一个候选状态$x'$，然后按照接受分布抽样决定是否接受状态$x'$（接受概率为$\alpha_{x,x'}$）。若接受，则时刻$t$转移到状态$x'$，否则仍停留在状态$x$。
• 对于建议分布的选择，主要有以下两种形式：

  1. 假设建议分布是对称的，即$q_{ij}=q_{ji},\forall i,j\in\mathcal{S}$。这种形式的建议分布也被称为Metropolis选择，是算法最初采用的形式。此时接受分布可简化为

     $$\alpha_{x,x'}=\min\left\{1,\frac{\pi_{x'}}{\pi_x}\right\}$$

     Metropolis选择有两种常用特例：

     • 第一种：建议分布取为条件分布$p(x'|x)$，其服从均值$x$，（协）方差为常数的（多元）正态分布；
     • 第二种：建议分布$q_{x,x'}=q(|x-x'|)$，这也被称为随机游走Metropolis。

     Metropolis选择的特点是当$x'$与$x$接近时，$q_{x,x'}$的概率值高，否则$q_{x,x'}$的概率值低。

  2. 假设建议分布与当前状态$x$无关，即$q_{x,x'}=q_{x'}$，也称为独立抽样。此时接受分布可写为

     $$\alpha_{x,x'}=\min\left\{1,\frac{\pi_{x'}}{q_{x'}}\cdot\frac{q_{x}}{\pi_{x}}\right\}$$

     独立抽样的实现更简单，但可能收敛速度慢，通常选择接近目标分布$\pi_x$的分布作为建议分布$q_x$。

• 当状态$x$为多元随机变量时，直接抽样往往比较困难，此时可以考虑对多元变量的每一分量的条件分布依次分别进行抽样，从而实现对整个多元变量的一次抽样，这种方法也被称作单分量Metropolis-Hastings算法。【具体略，可参见博客文章】

吉布斯抽样

在阐述吉布斯抽样之前，先介绍满条件分布：

• 设$\mathbf{x} = (x_1, x_2, \cdots, x_k)^\top$为$k$维随机变量，对应的马尔科夫链平稳分布$\pi_{\mathbf{x}} = \pi_{(x_1, x_2, \cdots, x_k)}$。如果条件概率分布$p(\mathbf{x}_I | \mathbf{x}_{-I})$中所有$k$个分量全部出现，其中 $$\mathbf{x}_I = \{x_i, i \in I\},\,\mathbf{x}_{-I} = \{x_i, i \notin I\},\quad I \subseteq K = \{1, 2, \cdots, k\}$$ 那么称这种条件概率分布为满条件分布。
  • 满条件分布有以下性质：对任意的$\mathbf{x} \in \mathcal{X}$和任意的$I \subseteq K$，有 $$p(\mathbf{x}_I | \mathbf{x}_{-I}) = \frac{\pi_{\mathbf{x}}}{\displaystyle\int \pi_{\mathbf{x}} \, \mathrm{d}\mathbf{x}_I} \propto \pi_{\mathbf{x}}$$ 因此对任意的$\mathbf{x}, \mathbf{x}' \in \mathcal{X}$和任意的$I \subseteq K$，有 $$\frac{p(\mathbf{x}'_I | \mathbf{x}'_{-I})}{p(\mathbf{x}_I | \mathbf{x}_{-I})} = \frac{\pi_{\mathbf{x}'}}{\pi_{\mathbf{x}}}$$ 这可以简化上述Metropolis-Hastings算法的运算。
• 吉布斯抽样的基本想法是：从联合概率分布定义满条件概率分布，依次对满条件概率分布进行抽样，得到样本的序列。
  • 具体而言，假设多元变量的联合概率分布为$p(\mathbf{x}) = p(x_1, x_2, \cdots, x_k)$。吉布斯抽样从一个初始样本$\mathbf{x}^{(0)} = (x_1^{(0)}, x_2^{(0)}, \cdots, x_k^{(0)})^\top$出发，不断进行迭代，每一次迭代得到联合分布的一个样本$\mathbf{x}^{(i)} = (x_1^{(i)}, x_2^{(i)}, \cdots, x_k^{(i)})^\top$。最终得到样本序列$\{\mathbf{x}^{(0)}, \mathbf{x}^{(1)}, \cdots, \mathbf{x}^{(n)}\}$。
  • 在每次迭代中，依次对$k$个分量中的一个分量进行随机抽样。如果在第$i$次迭代中，对第$j$个分量进行随机抽样，那么抽样的分布是满条件概率分布$p(x_{j}|\mathbf{x}_{-j}^{(i)})$，这里$\mathbf{x}_{-j}^{(i)}$表示第$i$次迭代中分量$j$以外的其他分量。
• 实际上，吉布斯抽样可以看作单分量Metropolis-Hastings算法的特殊情况。吉布斯抽样将建议分布取为满条件分布$q_{\mathbf{x},\mathbf{x}'}=p(x_j'|\mathbf{x}_{-j})$，而接受概率为 $$\begin{aligned} \alpha(\mathbf{x}, \mathbf{x}') &= \min\left\{1, \frac{\pi_{\mathbf{x}'} q_{\mathbf{x}', \mathbf{x}}}{\pi_{\mathbf{x}} q_{\mathbf{x}, \mathbf{x}'}}\right\}\\ &= \min\left\{1, \frac{p(x'_j | \mathbf{x}'_{-j}) p(x_j | \mathbf{x}'_{-j})}{p(x_j | \mathbf{x}_{-j}) p(x'_j | \mathbf{x}_{-j})}\right\} = 1 \end{aligned}$$ 其中利用了上述满条件分布的性质以及$\mathbf{x}_{-j}=\mathbf{x}'_{-j}$（每次只更新第$j$个分量）。所以其对应马尔科夫链的转移核就是满条件概率分布，吉布斯抽样对每次抽样的结果都接受，没有拒绝。

  此处默认$q_{\mathbf{x},\mathbf{x}'}$不为$0$，即马尔科夫链不可约。

• 下面给出吉布斯抽样的完整算法：
  1. 初始化。给出初样本$\mathbf{x}^{(0)} = (x_1^{(0)}, x_2^{(0)}, \cdots, x_k^{(0)})^\top$。
  2. 对$i$循环执行：
     设第$(i-1)$次迭代结束时的样本为$\mathbf{x}^{(i-1)} = (x_1^{(i-1)}, x_2^{(i-1)}, \cdots, x_k^{(i-1)})^\top$，则第$i$次迭代进行如下几步操作：
     1. 由条件分布$p(x_1 | x_2^{(i-1)}, x_3^{(i-1)}, \cdots, x_k^{(i-1)})$抽取$x_1^{(i)}$；
        $\vdots$
     2. 由条件分布$p(x_j | x_1^{(i)}, x_2^{(i)}, \cdots, x_{j-1}^{(i)}, x_{j+1}^{(i-1)}, \cdots, x_k^{(i-1)})$抽取$x_j^{(i)}$；
        $\vdots$
     3. 由条件分布$p(x_k | x_1^{(i)}, x_2^{(i)}, \cdots, x_{k-1}^{(i)})$抽取$x_k^{(i)}$；得到第$i$次迭代值$\mathbf{x}^{(i)} = (x_1^{(i)}, x_2^{(i)}, \cdots, x_k^{(i)})^\top$。
  3. 得到样本集合$\{\mathbf{x}^{(m+1)}, \mathbf{x}^{(m+2)}, \cdots, \mathbf{x}^{(n)}\}$。
  4. 计算样本均值 $$f_{mn} = \frac{1}{n-m} \sum_{i=m+1}^n f(\mathbf{x}^{(i)}).$$
• 吉布斯抽样适合于满条件概率分布容易抽样的情况，而单分量Metropolis-Hastings算法适合于满条件概率分布不容易抽样的情况，这时使用易于抽样的条件分布作为建议分布。
• 另一方面，对于一些概率图模型，利用条件独立性，可以将联合概率分布的吉布斯采样分解为条件概率分布的乘积的抽样。这样可以大幅减少抽样的计算复杂度。

---

ok，这样我们终于把无监督学习的内容整理完了！最后用两个表格简单进行总结：

---

顺利在开学前整完！敬请期待《深度学习》篇！