吉布斯采样

References：

• 浅谈「Gibbs采样」
• 走进贝叶斯统计（六）—— 吉布斯抽样 （Gibbs Sampling）
• 满条件分布
• 贝叶斯统计 | 第六章第二部分 二阶段吉布斯抽样
• 贝叶斯统计 | 第六章第三部分 多阶段吉布斯抽样

【概述】

当随机变量 $X$ 服从多维目标分布的情况下，对这个多维目标分布抽样，一种方法是 单分量 Metropolis Hasting 算法 中介绍的单分量 MH 算法

另一种则是吉布斯吉布斯采样（Gibbs Sampling），其基本思路是：从联合概率分布定义满条件概率分布，然后依次对满条件概率分布进行抽样，进而得到样本序列

单分量 MH 算法适合于满条件概率分布不容易抽样的情况，此时使用容易抽样的条件分布作为提议分布，而吉布斯采样适合于满条件概率容易抽样的情况

【满条件分布】

MCMC 的目标分布通常是多元联合概率分布 $p(\mathbf{x})=p(x_1,x_2,\cdots,x_k)$，其中 $\mathbf{x}=(x_1,x_2,\cdots,x_k)^T$ 是 $k$ 维随机变量

令

$$\begin{align*} \mathbf{x}_I=\{x_i,i\in I\} \\ \mathbf{x}_{-I}=\{x_i,i\notin I\} \\ I\subset K={1,2,\cdots,k} \end{align*}$$

若条件概率分布 $p(\mathbf{x}_I|\mathbf{x}_{-I})$ 中所有的 $k$ 个变量均出现，那么称该条件概率分布为满条件分布（Full Conditional Distribution）

满条件分布具备如下性质：

对任意的 $\mathbf{x},\mathbf{x}’\in\mathcal{X}$ 和 $I\subset K$，有：

$$p(\mathbf{x}_I|\mathbf{x}_{-I}) = \frac{p(\mathbf{x})}{\int p(\mathbf{x}) d\mathbf{x}_I} \propto p(\mathbf{x})$$

且满足：

$$\frac{p(\mathbf{x}'_{I} | \mathbf{x}_{-I}' )}{p(\mathbf{x}_{I}|\mathbf{x}_{-I})} = \frac{p(\mathbf{x}')}{p(\mathbf{x})}$$

吉布斯采样正是通过满条件分布来对每一个维度的分量进行抽样

【基本思路】

假设多元随机变量的联合概率分布为 $p(\mathbf{x})=p(x_1,x_2,\cdots,x_k)$，吉布斯采样从一个初始样本 $\mathbf{x}^{(0)}=(x_1^{(0)},x_2^{(0)},\cdots,x_k^{(0)})^T$ 出发，不断进行迭代，每一次迭代得到联合分布的一个样本 $\mathbf{x}^{(i)}=(x_1^{(i)},x_2^{(i)},\cdots,x_k^{(i)})^T$，最终得到样本序列 $\{\mathbf{x}^{(0)},\mathbf{x}^{(2)},\cdots,\mathbf{x}^{(n)}\}$

在每次迭代中，依次对 $k$ 个随机变量中的一个变量进行随机抽样，若在第 $i$ 次迭代中，对第 $j$ 个变量进行随机抽样，那么抽样的分布是满条件概率分布 $p(x_j|x_{-j}^{(i)})$，其中 $x_{-j}^{(i)}$ 代表第 $i$ 次迭代中，分量 $j$ 以外的其他变量

设在第 $i-1$ 步得到样本 $(x_1^{(i-1)},x_2^{(i-1)},\cdots,x_k^{(i-1)})^T$，在第 $i$ 步，按照如下步骤进行随机抽样：

Step 1：对第一个变量按以下满条件概率分布随机抽样得到 $x_1^{(i)}$

$$p(x_1|x_2^{(t-1)},\cdots,x_k^{(t-1)})$$

Step 2：依次对 $j$ 个变量按以下满条件概率分布随机抽样得到 $x_j^{(i)}$

$$p(x_j|x_1^{(i)},\cdots,x_{j-1}^{(i)},x_{j-1}^{(i-1)},\cdots,x_k^{(i-1)}),\quad j=2,\cdots,k-1$$

Step 3：对第 $k$ 个变量按以下满条件概率分布随机抽样得到 $x_k^{(i)}$

$$p(x_k|x_1^{(i)},\cdots,x_{k-1}^{(i)})$$

故可得到整体样本 $\mathbf{x}^{(i)}=(x_1^{(i)},x_2^{(i)},\cdots,x_k^{(i)})^T$

【与 MH 算法的联系】

吉布斯采样是单分量 MH 算法的特殊情况

令转移核为满条件概率分布 $p(x’_j|x_{-j})\neq0$，即：

$$p(x,x')=p(x'_j|x_{-j})$$

并令提议分布是当前变量 $x_j$ 的满条件概率分布，即：

$$q(x,x')=p(x_j'|x_{-j})$$

由于：

$$p(x_{-j})=p(x'_{-j}) \quad p(\cdot|x_{-j})=p(\cdot|x'_{-j})$$

故对接受概率，有：

$$\begin{align*} \alpha(x,x') &= \min \{1, \frac{p(x')q(x',x)}{p(x)q(x,x')} \} \\ &= \min \{1, \frac{p(x')p(x_{j}|x_{-j}')}{p(x)p(x_j'|x_{-j})} \} \\ &= \min \{1, \frac{p(x'_{-j})p(x'_j|x_{-j})p(x_j|x'_{-j})}{p(x_{-j})p(x_{j}|x_{-j})p(x'_{j}|x_{-j})} \} \\ &= 1 \end{align*}$$

也就是说，依次按照单变量的满条件概率分布 $p(x’_j|x_{-j})$ 进行随机抽样，就能实现单分量 MH 算法

此外，需要注意的是，吉布斯采样对每次抽样的结果都接受，这一点与 MH 算法不同

【算法流程】

综上所述，吉布斯采样的算法流程如下

输入：目标概率分布的概率密度函数 $p(\mathbf{x})$，函数 $f(\mathbf{x})$

输出：$p(\mathbf{x})$ 的随机样本 $\{ \mathbf{x}^{(m+1)},\mathbf{x}^{(m+2)},\cdots,\mathbf{x}^{(n)} \}$，函数样本均值 $\hat{f}_{mn}$

参数：收敛步数 $m$，迭代步数 $n$

算法步骤：

Step 1：算法初始化，随机取得或人工定义第一个样本 $\mathbf{x}^{(0)}=(x_1^{(0)},x_2^{(0)},\cdots,x_k^{(0)})^T$

Step 2：对 $i=1,2,\cdots,n$ 循环执行

设第 $i-1$ 次迭代结束时的样本为 $\mathbf{x}^{(i-1)}=(x_1^{(i-1)},x_2^{(i-1)},\cdots,x_k^{(i-1)})^T$，则第 $i$ 次迭代进行如下几步操作

1）由满条件分布 $p(x_1|x_2^{(i-1)},x_3^{(i-1)},\cdots,x_k^{(i-1)})$ 抽取 $x_1^{(i)}$

2）由满条件分布 $p(x_2|x_1^{(i)},x_3^{(i-1)},\cdots,x_k^{(i-1)})$ 抽取 $x_2^{(i)}$

$$\vdots$$

j）由满条件分布 $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$$

k）由满条件分布 $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)})^T$

Step 3：从样本集合 $\{ \mathbf{x}^{(0)},\mathbf{x}^{(1)},\cdots,\mathbf{x}^{(k)} \}$ 中得到样本集合 $\{ \mathbf{x}^{(m+1)},\mathbf{x}^{(m+2)},\cdots,\mathbf{x}^{(n)} \}$，并计算函数样本均值 $\hat{f}_{mn}$

$$\hat{f}_{mn}=\frac{1}{n-m}\sum_{i=m+1}^{n} f(\mathbf{x}^{(i)})$$

【示例】

使用吉布斯采样，从以下二元正态分布中抽取随机样本

$$\mathbf{x}=(x_1,x_2)^T \sim p(x_1,x_2)=N(0,\Sigma), \quad \Sigma=\begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix}$$

对于 $p(x_1,x_2)\sim N(0,\Sigma)$，条件概率分布为一元正态分布：

$$\begin{align*} p(x_1|x_2) = N(\rho x_2,(1-\rho^2)) \\ p(x_2|x_1) = N(\rho x_1,(1-\rho^2)) \\ \end{align*}$$

假设初始样本 $\mathbf{x}^{(0)}=(x_1^{(0)},x_2^{(0)})$，通过吉布斯采样，可以得到如下图所示的样本序列：

得到的样本集合 $\{ \mathbf{x}^{(m+1)},\mathbf{x}^{(m+2)},\cdots,\mathbf{x}^{(n)} \}$ 就是二元正态分布得多随机抽样

下图展示了吉布斯采样的迭代过程

单分量 MH 算法在抽样时，在样本间移动时可能会在某些点上停留（抽样被拒绝），而吉布斯采样在抽样时，会在样本点间持续移动