什么是吉布斯采样（Gibbs Sampling）

---

1 蒙特卡洛方法

介绍吉布斯采样前，我们先看一下蒙特卡洛方法。

1.1 蒙特卡洛方法的作用

有很多函数我们无法直接得到他的积分值，但我们可以利用蒙特卡洛方法来进行估计。
比如下面的积分：
$${\int }_a^{b}f(x)dx$$
我们采样一个点作为函数$f(x)$的均值对积分进行估计：
$${\int }_a^{b}f(x)dx\approx f(x_0)(b-a)$$
只采样一个点，误差比较大，我们可等间隔的采样n个点：
$${\int }_a^{b}f(x)dx\approx \frac{b-a}{n}\sum _{i=0}^{n-1}f(x_i)$$

1.2 非均匀分布采样

上述我们是等间隔采样，其实隐含着假设：x在(a,b)区间是均匀分布的。如果x在(a,b)区间不是均匀分布，而是按照概率p(x)分布的呢？此时利用采样点进行积分估计的结果如下所示：
$${\int }_a^{b}f(x)dx={\int }_a^b\frac{f(x)}{p(x)}p(x)dx\approx \frac{1}{n}\sum _{i=0}^{n-1}\frac{f(x_i)}{p(x_i)}$$
相当于我按照概率p(x)采样n个点，然后根据这n个点对应的$f(x_i)$和$p(x_i)$值来进行估计。均匀分布相当于一种特例。如下：
$$p(x)=\frac{1}{b-a}$$
$$\begin{gathered} {\int }_a^{b}f(x)dx\approx \frac{1}{n}\sum _{i=0}^{n-1}\frac{f(x_i)}{p(x_i)} \\ =\frac{1}{n}\sum _{i=0}^{n-1}\frac{f(x_i)}{1/(b-a)} \\ =\frac{b-a}{n}\sum _{i=0}^{n-1}f(x_i) \end{gathered}$$

1.3 分布p(x)不好采样怎么办？

假设我们知道概率分布p(x)，但是这个概率分布不常见，不知道怎么采样。我们可以借助容易采样的分布q(x)来对p(x)进行采样。方法如下：

1. 首先选择一个k,使得$p(x)⩽kq(x)$；
2. 然后我们依照q(x)分布采样一个$x_0$；
3. 然后我们在均匀分布[0,kq(x_0)]区间随机选择一个点$z_0$；
4. 若$z_0⩽p(x_0)$，则接受这个采样点$x_0$，否则，就丢弃这个采样点；
5. 重复上述步骤即可得到符合p(x)分布的样本点。

2 什么是吉布斯采样

如果一个分布，既无法直接采样，也无法借助上文中的方法进行采样呢？如下：

1. 不知道分布$p(x,y)$，只知道其条件分布$p(x|y),p(y|x)$，无法利用上面的方法得到符合$p(x,y)$分布的样本点
2. 高维的概率分布$p(x_1,x_2,...,x_n)$，找不到分布q满足$p⩽kq$。

此时，我们就要引入马尔可夫链，借助马尔可夫链的一些特性来解决采样的问题。

2.1 马尔可夫链

2.1.1 什么是马尔可夫链呢？

简单来说就是一个状态转移过程，且当前状态只跟前一个状态有关。以股市为例，假设有牛市，熊市，横盘三种状态，牛市有0.3的概率保持牛市，有0.2的概率转为熊市，有0.5的概率横盘，这种转移过程就是马尔可夫链。
转移过程会有一个初始状态$\pi (x_o)$，会有一个状态空间，还有一个状态转移矩阵P。
转移过程中，若出现$\pi (x_{n-1})P=\pi (x_n)=\pi (x_{n-1})$ ，即$\pi P=\pi$，则称概率分布$\pi$是状态转移矩阵P的平稳分布。也称其为平稳的马尔可夫链。
马尔可夫链还有一些比较好的性质：
性质：给定初始状态$\pi (x_0)$，给定转移矩阵P，n轮之后，$\pi$会收敛。且稳定后的$\pi$与初始状态$\pi (x_0)$无关。
性质：对于给定的转移矩阵$P$,$P^n$在n大于一定值的时候是确定的。

2.1.2 为什么我们要引入马尔可夫链？

其实马尔可夫链真正对我们有用的是这个公式：$\pi P=\pi$
其中$\pi$是我们想要采样的分布，P是一个状态转移矩阵。
我们只要找到这样的P，就可以把对分布$\pi (x)$进行采样，转化为一个转移过程，因为转移过程中生成的样本点，就是我们想要的符合$\pi (x)$分布的样本点。

2.1.3 对给定的分布$\pi$，怎么找到对应的P，使得其为平稳马尔可夫过程

有一个细致平稳条件：
$$\pi (i)P(i,j)=\pi (j)P(j,i)$$

由细致平稳条件，我们就能得到$\pi P=\pi$，推导如下：
$$\sum _i\pi (i)P(i,j)=\sum _i\pi (j)P(j,i)=\pi (j)\sum _{i}P(j,i)=\pi (j)$$
上式就是$\pi P=\pi$ 元素展开形式。

所以我们只要找到满足细致平稳条件的P，就能利用马尔可夫链对$\pi$采样了。

2.2 MCMC采样

假设$\pi (x)$为我们期望进行采样的目标分布，$Q$为状态转移矩阵，其不一定满足细致平稳条件$\pi (i)Q(i,j)=\pi (j)Q(j,i)$.

为了使其满足细致平稳条件，我们引入一个$\alpha (i,j)$,使得下式成立：
$$\pi (i)Q(i,j)\alpha (i,j)=\pi (j)Q(j,i)\alpha (j,i)$$
此时对应的 ，状态转移矩阵$P(i,j)=Q(i,j)\alpha (i,j)$。

怎么才能找到这样的$\alpha (i,j)$呢？其实只要$\alpha (i,j)$满足下面条件即可：
$$\begin{gathered} \alpha (i,j)=\pi (j)Q(j,i) \\ \alpha (j,i)=\pi (i)Q(i,j) \end{gathered}$$
上面两个式子是等价的，我们只需直接令$\alpha (i,j)=\pi (j)Q(j,i)$就行了。

其中，$\alpha (i,j)$我们称之为接受率。有点类似蒙特卡洛采样那里的接受-拒绝方法。

MCMC采样过程如下：

1. 已知状态转移矩阵Q，分布$\pi (x)$。设状态转移次数阈值为$n_1$，需要样本个数为$n_2$。
2. 从任意指定初始状态$x_0$,令 $t=0$.
3. while $t<n_1+n_2$
   a. 从条件分布$Q(x|x_t)$中采样得到样本$x_{\ast }$
   b. 从[0,1]均匀分布采样得到u
   c. 若$u\le \alpha (x_t,x_{\ast })=\pi (x_{\ast })Q(x_{\ast },x_t)$, 则接受转移，即$x_{t+1}=x_{\ast }$，t++
   d. 否则，不接受转移。t不变。
4. 最后得到的样本集$x_{n_1},x_{n_1+1},...,x_{n_1+n_2-1}$就是我们需要的采样结果。

2.3 M-H采样

M-H采样是Metropolis Hastings采样的简称。

MCMC已经解决了我们的采样问题，但是其样本接受率可能比较低，就是$\alpha (x_t,x_{\ast })$可能为0.1，0.001，这样采样次数就会比较多，浪费了很多时间。怎样提高接受率呢？
观察细致平稳条件
$$\pi (i)Q(i,j)\alpha (i,j)=\pi (j)Q(j,i)\alpha (j,i),其中\alpha (i,j)=\pi (j)Q(j,i)$$
我们知道，在两边同时乘以一个数字，或者除以一个非0数字，细致平稳条件也是满足的。所以对上式两边同乘 min ⁡ { α ( j , i ) , α ( i , j ) } \min \{\alpha(j,i),\alpha(i,j)\} min{α(j,i),α(i,j)}，可得下式：
π ( i ) Q ( i , j ) α ( i , j ) ∗ min ⁡ { α ( j , i ) , α ( i , j ) } = π ( j ) Q ( j , i ) α ( j , i ) ∗ min ⁡ { α ( j , i ) , α ( i , j ) } \pi (i)Q(i,j)\alpha(i,j)*\min \{\alpha(j,i),\alpha(i,j)\}=\pi (j)Q(j,i)\alpha(j,i) *\min \{\alpha(j,i),\alpha(i,j)\} π(i)Q(i,j)α(i,j)∗min{α(j,i),α(i,j)}=π(j)Q(j,i)α(j,i)∗min{α(j,i),α(i,j)}
然后两边再同时除以$\alpha (j,i)\ast \alpha (i,j)$,得：
π ( i ) Q ( i , j ) ∗ min ⁡ { 1 , α ( i , j ) α ( j , i ) } = π ( j ) Q ( j , i ) ∗ min ⁡ { α ( j , i ) α ( i , j ) , 1 } \pi (i)Q(i,j)*\min \{1,\frac{\alpha(i,j)}{\alpha(j,i)}\}=\pi (j)Q(j,i)*\min \{\frac{\alpha(j,i)}{\alpha(i,j)},1\} π(i)Q(i,j)∗min{1,α(j,i)α(i,j)​}=π(j)Q(j,i)∗min{α(i,j)α(j,i)​,1}
所以我们得到：

α n e w ( i , j ) = min ⁡ { α ( i , j ) α ( j , i ) , 1 } = min ⁡ { π ( j ) Q ( j , i ) π ( i ) Q ( i , j ) , 1 } \alpha_{new}(i,j)=\min \{\frac{\alpha(i,j)}{\alpha(j,i)},1\}=\min \{ \frac{\pi (j)Q(j,i)}{\pi (i)Q(i,j)},1\} αnew​(i,j)=min{α(j,i)α(i,j)​,1}=min{π(i)Q(i,j)π(j)Q(j,i)​,1}

如果转移矩阵Q是对称的，上式可简化为：
α n e w ( i , j ) = min ⁡ { π ( j ) π ( i ) , 1 } \alpha_{new}(i,j)=\min \{ \frac{\pi (j)}{\pi (i)},1\} αnew​(i,j)=min{π(i)π(j)​,1}

M-H采样是对MCMC方法的改进，也属于MCMC采样。

2.4 吉布斯采样（Gibbs）

但是M-H采样也有一些问题，比如

• 当特征比较多时，$\pi (j)Q(j,i)/\pi (i)Q(i,j)$的计算比较耗时。
• 虽然${\alpha }_{new}(i,j)$的接受率比较大，但还是小于1的，还是有很多辛苦计算的采样点被拒绝掉。
• 此外，很多时候我们不知道其联合分布，只知道其条件分布，那么这种情况我们就无法对联合分布进行采样。

吉布斯采样主要就是为了解决了上述问题。

• 吉布斯采样也是一种MCMC采样方法，可以看作是M-H算法的一个特例。
• 吉布斯采样适用于联合分布未明确知道或难以直接抽样，但每个变量的条件分布是已知的，并且很容易（或者至少更容易）从中抽样的情况。
• 吉布斯采样之所以被看作是M-H算法的一个特例，是因为它总是以1的概率接受抽样出的值。
• 吉布斯采样要求数据至少是两维的。M-H采样无此要求。

2.4.1 吉布斯采样原理

在MCMC采样中，我们要寻找一个满足细致平稳条件$\pi (i)P(i,j)=\pi (j)P(j,i)$的状态转移矩阵P。

下面，我们利用条件概率的变换，来寻找这样的P。首先我们以见到那的二维情况维例：

2.4.1.1 二维情况

假设我们有状态$A=(x_1^1,x_2^1)$, $B=(x_1^1,x_2^2)$, $C=(x_1^2,x_2^1)$, 则有：
$$\begin{gathered} \pi (A)\pi (x_2^2|x_1^1)=\pi (x_1^1,x_2^1)\pi (x_2^2|x_1^1) \\ =\pi (x_1^1)\pi (x_2^1|x_1^1)\pi (x_2^2|x_1^1) \\ =\pi (x_1^1)\pi (x_2^2|x_1^1)\pi (x_2^1|x_1^1) \\ =\pi (x_1^1,x_2^2)\pi (x_2^1|x_1^1) \\ =\pi (B)\pi (x_2^1|x_1^1) \end{gathered}$$
$$\begin{gathered} \pi (A)\pi (x_1^2|x_2^1)=\pi (x_1^1,x_2^1)\pi (x_1^2|x_2^1) \\ =\pi (x_2^1)\pi (x_1^1|x_2^1)\pi (x_1^2|x_2^1) \\ =\pi (x_2^1)\pi (x_1^2|x_2^1)\pi (x_1^1|x_2^1) \\ =\pi (x_1^2,x_2^1)\pi (x_1^1|x_2^1) \\ =\pi (C)\pi (x_1^1|x_2^1) \end{gathered}$$

如果我们令状态转移概率矩阵P为下式：

$P(A\to B)=\pi (x_2^B|x_1^1)，若x_1^A=x_1^B=x_1^1$
$P(A\to C)=\pi (x_1^C|x_2^1)，若x_2^A=x_2^C=x_2^1$
$P(A\to D)=0，其他情况$

那么对任意两点E,F, 将满足细致平稳条件。
$$\pi (E)P(E\to F)=\pi (F)P(F\to E)$$

证明如下：
因为满足
$x_1^A=x_1^B=x_1^1$，
所以
$$\begin{gathered} P(A\to B)=\pi (x_2^B|x_1^1)=\pi (x_2^2|x_1^1) \\ P(B\to A)=\pi (x_2^A|x_1^1)=\pi (x_2^1|x_1^1) \end{gathered}$$
所以有：
$\pi (A)P(A\to B)=\pi (A)\pi (x_2^2|x_1^1)=\pi (B)\pi (x_2^1|x_1^1)=\pi (B)P(B\to A)$

因为满足
$x_2^A=x_2^C=x_2^1$
所以
$P(A\to C)=\pi (x_1^C|x_2^1)=\pi (x_1^2|x_2^1)$
$P(C\to A)=\pi (x_1^A|x_2^1)=\pi (x_1^1|x_2^1)$
所以有：
$\pi (A)P(A\to C)=\pi (A)\pi (x_1^2|x_2^1)=\pi (C)\pi (x_1^1|x_2^1)=\pi (C)\pi (C\to A)$

因为既不满足
$x_1^A=x_1^D=x_1^1$
也不满足
$x_2^A=x_2^D=x_2^1$
所以
$P(A\to D)=0$
$P(D\to A)=0$
所以有：
$\pi (A)P(A\to D)=0=\pi (D)\pi (D\to A)$

综上，得证任意两点E,F, 有
$\pi (E)P(E\to F)=\pi (F)P(F\to E)$

2.4.1.2 高维情况

假设假设$A=(x_1^1,x_2^1,...,x_n^1)$,$B_n=(x_1^1,x_2^1,...x_{n-1}^1,x_n^2)$,我们有：
$$\begin{gathered} \pi (A)\pi (x_n^2|x_1^1,x_2^1,...,x_{n-1}^1)=\pi (x_1^1,x_2^1,...,x_n^1)\pi (x_n^2|x_1^1,x_2^1,...,x_{n-1}^1) \\ =\pi (x_1^1,x_2^1,...,x_{n-1}^1)\pi (x_n^1|x_1^1,x_2^1,...,x_{n-1}^1)\pi (x_n^2|x_1^1,x_2^1,...,x_{n-1}^1) \\ =\pi (x_1^1,x_2^1,...,x_{n-1}^1)\pi (x_n^2|x_1^1,x_2^1,...,x_{n-1}^1)\pi (x_n^1|x_1^1,x_2^1,...,x_{n-1}^1) \\ =\pi (x_1^1,x_2^1,...x_{n-1}^1,x_n^2)\pi (x_n^1|x_1^1,x_2^1,...,x_{n-1}^1) \\ =\pi (B)\pi (x_n^1|x_1^1,x_2^1,...,x_{n-1}^1) \end{gathered}$$

因此，如果我们令状态转移概率矩阵为下式：

• $P(A\to B_n)=\pi (x_n^{B_n}|x_1^1,x_2^1,...,x_{n-1}^1)，若x_i^A=x_i^{B_n}=x_i^1,1\le i\le n且i\ne n.$
• $P(A\to B_{n-1})=\pi (x_{n-1}^{B_{n-1}}|x_1^1,x_2^1,...,x_{n-2}^1,x_n^1)，若x_i^A=x_i^{B_{n-1}}=x_i^1,1\le i\le n且i\ne n-1.$
• …
• $P(A\to B_j)=\pi (x_j^{B_j}|x_1^1,x_2^1,...,x_{j-1}^1,x_{j+1}^1,...,x_n^1)，若x_i^A=x_i^{B_j}=x_i^1,1\le i\le n且i\ne j.$
• …
• $P(A\to B_1)=\pi (x_1^{B_1}|x_2^1,x_3^1,...,x_n^1)，若x_i^A=x_i^{B_1}=x_i^1,1\le i\le n且i\ne 1.$
• $P(A\to D)=0，其他情况.$

那么对任意两点E，F 将满足细致平稳条件,$\pi (E)P(E\to F)=\pi (F)P(F\to E)$

证明过程同二维情况。

---

综上，状态转移矩阵P就变成了完全条件概率：
$A=(x_1^1,x_2^1,...,x_n^1)$,
$B_j=(x_1^1,x_2^1,...x_{j-1}^1,x_j^2,x_{j+1}^1,...,x_n^1)$
$P(A\to B_j)=\pi (x_j^{B_j}|x_1^1,x_2^1,...,x_{j-1}^1,x_{j+1}^1,...,x_n^1)$
此时，我们的采样过程就变成了按轴转移的过程了。想更新样本点的第j维，使用第j维的完全条件概率分布$\pi (x_j^{new}|x_1^{old},x_2^{old},...,x_{j-1}^{old},x_{j+1}^{old},...,x_n^{old})$,即可采样得到第j维新的值$x_j^{new}$。
所有维度的值都更新一遍，就能得到一个新的样本点$(x_1^{new},x_2^{new},...,x_n^{new})$

2.4.2 吉布斯采样过程

具体采样过程如下：

1. 随机初始化 { x i 0 : i = 1 , 2 , . . . , M } \{x_i^0:i=1,2,...,M\} {xi0​:i=1,2,...,M}
2. for $\tau =0,1,2,...,n_1+n_2-1$，其中$n_1$是采样阈值，$n_2$是需要的样本个数。
   采样$x_1^{\tau +1}\sim p(x_1|x_2^{\tau },x_3^{\tau },...,x_M^{\tau })$
   采样$x_2^{\tau +1}\sim p(x_2|x_1^{\tau +1},x_3^{\tau },...,x_M^{\tau })$
   …
   采样$x_j^{\tau +1}\sim p(x_j|x_1^{\tau +1},x_2^{\tau +1},...,x_{j-1}^{\tau +1},x_{j+1}^{\tau }...,x_M^{\tau })$
   …
   采样$x_M^{\tau +1}\sim p(x_M|x_1^{\tau +1},x_2^{\tau +1},...,x_{M-1}^{\tau +1})$
3. 以此类推，得到采样点$(x_1^{n_1},x_2^{n_1},...,x_M^{n_1}),(x_1^{n_1+1},x_2^{n_1+1},...,x_M^{n_1+1}),...,(x_1^{n_1+n_2-1},x_2^{n_1+n_2-1},...,x_M^{n_1+n_2-1})$。

上面过程是依次轮换坐标轴进行采样。
实际上，按照随机顺序变换坐标轴也是可以的。

参考资料

该文主要参考下面四篇文章，并加上一些自己的理解和推理细节，感谢作者，写的真好：
MCMC(一)蒙特卡罗方法
MCMC(二)马尔科夫链
MCMC(三)MCMC采样和M-H采样
MCMC(四)Gibbs采样