在定义了一个贝叶斯模型后，可以从结果的后验分布中抽样来对参数进行推断。除非是简单的模型，这通常是个复杂的计算过程。本节介绍用 马尔可夫链蒙特卡洛（Markov chain Monte Carlo，MCMC） 的方法进行后验分布抽样。并可将吉布斯采样（Gibbs sampling），一种 MCMC 流程，与最大期望（EM）算法相联系：两者主要的不同是一个从条件分布中抽样而另一个在条件分布上求最大化解。
考虑这样一个抽象问题。有若干随机变量 $U_1$，$U_2$，……，$U_K$，我们想要从它们的联合分布进行抽样。假设这个过程并不简单，但模拟条件分布 $\operatorname{Pr}(U_j|U_1,U_2,\dots,U_{j-1},U_{j+1},\dots,U_K)$，$j=1,2,\dots,K$要相对容易。吉布斯采样过程就是在模拟每个条件分布，并且当过程结果稳定后就得到了一个目标联合分布的采样。算法 8.3 描述了这个过程。
算法 8.3：吉布斯采样（Gibbs sampling）
- 选取某些初始值 $U_k^{(0)}$，$k=1,2,\dots,K$。
- 对 $t=1,2,\dots$ 循环进行：
依次为 $k=1,2,\dots,K$，根据条件分布 $\operatorname{Pr}(U_k^{(t)}|U_1^{(t)},\dots,U_{k-1}^{(t)},U_{k+1}^{(t)},\dots,U_{K}^{(t)})$ 生成 $U_k^{(t)}$。 - 重复步骤 2 直到 $(U_1^{(t)},U_2^{(t)},\dots,U_K^{(t)})$ 的联合分布不再改变。
在一些正则条件下，可证明这个过程最终趋于稳定，并且最终的随机变量确实是 $(U_1,U_2,\dots,U_K)$ 联合分布的一次抽样。尽管样本 $(U_1^{(t)},U_2^{(t)},\dots,U_K^{(t)})$ 明显会根据 $t$ 变化。更正式解释是吉布斯采样过程产生了一个以真实联合分布为稳定分布的马尔可夫链，这也是“马尔可夫链蒙特卡洛”名字的含义。因为后续步骤不会改变 $U_k$ 的边际分布，在这个过程中真实联合分布自然与稳定分布相同1。
注意在过程中并不需要条件分布的显式表达式，而只需要能够从中进行抽样。在过程到达稳态后，任意变量子集的边际密度函数可近似为利用样本取值的密度函数估计。然而，如果已知条件密度 $\operatorname{Pr}(U_k|U_l,l\ne k)$ 的显式表达式，以 $U_k$ 的边际密度函数为例，则一个更好的估计为（练习 8.3）：
$$\widehat{\operatorname{Pr}}_{U_k}(u) = \frac{1}{M-m+1} \sum_{t=m}^M \operatorname{Pr}(u|U_\ell^{(t)}, \ell \ne k) \tag{8.50}$$其中对序列的后 $M-m+1$ 个元素取平均，这是考虑到过程在达到稳态之前的初始“预热”阶段。
再回到贝叶斯推断，我们的目的是从给定数据 $\mathbf{Z}$ 所得到的参数联合后验分布中抽取一个样本。如果从每个参数的条件于给定其他参数和数据 $\mathbf{Z}$ 上的分布中更容易抽样，那么就可以利用吉布斯采样。下面描述在高斯混合分布问题中的具体过程。
从后验分布的吉布斯采样与指数分布族模型的最大期望（EM）算法有紧密的联系。关键在于将最大期望算法中的隐变量 $\mathbf{Z}^m$ 看做是吉布斯采样的另一个参数。以高斯混合问题为例，现在的参数为 $(\theta,\mathbf{Z}^m)$。简单起见，将方差 $\sigma_1^2$ 和 $\sigma_2^2$ 以及混合比重 $\pi$ 固定在他们的最大似然值上，因此 $\theta$ 中的未知参数只有均值 $\mu_1$ 和 $\mu_2$。算法 8.4 给除了混合模型问题中的吉布斯采样。可见步骤 2.1 和 2.2 与最大期望算法中的 E 和 M 步骤想对应，只是一个进行抽样一个进行最大化。在步骤 2.1 中，吉布斯采样没有计算最大似然责任值 $\gamma_i=\operatorname{E}(\Delta_i|\theta,\mathbf{Z})$，而是从概率 $\operatorname{Pr}(\Delta_i|\theta,\mathbf{Z})$ 中模拟隐变量 $\Delta_i$。在步骤 2.2 中，并没有对后验概率 $\operatorname{Pr}(\mu_1,\mu_2,\mathbf{\Delta}|\mathbf{Z})$ 最大化，而是从条件概率 $\operatorname{Pr}(\mu_1,\mu_2|\mathbf{\Delta},\mathbf{Z})$ 中模拟数据。
算法 8.4：混合模型的吉布斯采样
- 选取某些初始值 $\theta^{(0)}=(\mu_1^{(0)},\mu_2^{(0)})$
- 对 $t=1,2,\dots$ 循环进行：
- 基于概率 $\operatorname{Pr}(\Delta_i=1)=\hat{\gamma}_i(\theta^{(t)})$ 生成 $\Delta_i\in\{0,1\}$，其中的概率来自表达式 8.42： $$\hat{\gamma}_i = \frac{\hat{\pi}\phi_{\hat{\theta}_2}(y_i)} {(1-\hat{\pi})\phi_{\hat{\theta}_1}(y_i) + \hat{\pi}\phi_{\hat{\theta}_2}(y_i)}, i=1,2,\dots,N \tag{8.42}$$
- 计算 $$\begin{align} \hat{\mu}_1 &= \frac{\sum_{i=1}^N (1-\Delta_i^{(t)}) \cdot y_i} {\sum_{i=1}^N (1-\Delta_i^{(t)})} \\ \hat{\mu}_2 &= \frac{\sum_{i=1}^N \Delta_i^{(t)} \cdot y_i} {\sum_{i=1}^N \Delta_i^{(t)}} \end{align}$$ 并基于分布生成 $\mu_1\sim\mathcal{N}(\hat{\mu}_1,\hat{\sigma}_1^2)$ 和 $\mu_2\sim\mathcal{N}(\hat{\mu}_2,\hat{\sigma}_2^2)$。
- 重复步骤 2 直到 $(\mathbf{\Delta}^{(t)},\mu_1^{(t)},\mu_2^{(t)})$ 的联合分布不再改变。
图 8.8 展示了重复了 200 次的吉布斯采样，左图的上下分别为均值参数 $\mu_1$ 和 $\mu_2$，右图为观测样本中类别 2 的比例 $\sum_i\Delta_i/N$。每个图中的水平虚线分别为最大似然估计值 $\hat{\mu}_1$，$\hat{\mu}_2$ 和 $\sum_i\hat{\gamma}_i/N$。可见抽样值很快达到稳定，并且均匀地分布在最大似然值的周围。
上述为了阐明吉布斯采样和最大期望算法简化了混合模型。在更现实的条件下，可能要假设方差 $\sigma_1^2$ 和 $\sigma_2^2$ 以及混合比例 $\pi$ 的先验分布，并在吉布斯采样中添加单独的从条件于其他参数下它们的后验分布抽样的步骤。需要选择合适的先验分布，否则有可能得到一个退化的后验分布，即全部的混合权重都集中在一个成分上。
吉布斯采样是众多近期提出的后验分布抽样方法中的一个。它使用了条件于其他参数下每个参数的条件分布进行抽样，可用于模型结构恰好容易进行条件分布抽样的问题中。有不要求这种结构的其他方法，例如 **梅特罗波利斯－黑斯廷斯（Metropolis-Hastings）**算法。这些以及其他计算的贝叶斯方法用于复杂学习方法中，比如高斯过程模型和神经网络。更多细节可参阅本章最后的参考文献说明。
本节练习
练习 8.3
Justify the estimate (8.50), using the relationship
$$\operatorname{Pr}(A) = \int \operatorname{Pr}(A|B) d\operatorname{Pr}(B)$$-
不理解这句话：It is not surprising that the true joint distribution is stationary under this process, as the successive steps leave the marginal distributions of the $U_k$’s unchanged. ↩︎