引言

上一节介绍了将马尔可夫链与蒙特卡洛方法相结合的算法——MH采样算法(Metropolis Hastings)，本节将介绍吉布斯采样算法(Gibbs Sampling)。

回顾：MH采样算法

基于马尔可夫链的采样方式

首先，MH算法是基于马尔可夫链的采样方式：

通过构建合适的状态转移矩阵A∗\mathcal A^*A∗，使得连续状态下的随机变量Xt,Xt+1\mathcal X_t,\mathcal X_{t+1}Xt,Xt+1对应的概率分布满足细致平衡原则(Detail Balance)：
π(Xt=xi)⋅P(Xt+1=xj∣Xt=xi)=π(Xt=xj)⋅P(Xt+1=xi∣Xt=xj)\pi(\mathcal X_{t} =x_i) \cdot P(\mathcal X_{t+1} =x_j \mid \mathcal X_t = x_i) = \pi(\mathcal X_{t} = x_j) \cdot P(\mathcal X_{t+1} =x_i \mid \mathcal X_t =x_j)π(Xt=xi)⋅P(Xt+1=xj∣Xt=xi)=π(Xt=xj)⋅P(Xt+1=xi∣Xt=xj)
其中，π(Xt),π(Xt+1)\pi(\mathcal X_t),\pi(\mathcal X_{t+1})π(Xt),π(Xt+1)表示不同时刻对应随机变量Xt,Xt+1\mathcal X_t,\mathcal X_{t+1}Xt,Xt+1的概率分布；P(Xt+1∣Xt)P(\mathcal X_{t+1} \mid \mathcal X_t)P(Xt+1∣Xt)表示ttt时刻到t+1t+1t+1时刻的状态转移概率。
这种状态转移矩阵A∗\mathcal A^*A∗的构建方式，使得马尔可夫链{XT}\{\mathcal X_{T}\}{XT}共用同一个概率分布，即平稳分布；
已知上一状态随机变量Xk−1\mathcal X_{k-1}Xk−1选择具体离散结果xk−1x_{k-1}xk−1的条件下，当前状态随机变量Xk\mathcal X_{k}Xk的采样结果可表示如下：
xk∼P(Xk∣Xk−1=xk−1)x_k \sim P(\mathcal X_{k} \mid \mathcal X_{k-1} = x_{k-1})xk∼P(Xk∣Xk−1=xk−1)
由于平稳分布的原因，使得 各时间状态产生的样本点{x0,x1,⋯,xk−1,xk,⋯}\{x_0,x_1,\cdots,x_{k-1},x_{k},\cdots \}{x0,x1,⋯,xk−1,xk,⋯}均属于平稳分布的样本点，从而使用蒙特卡洛方法来求解该平稳分布的期望信息：
NNN表示采集样本点数量。
EX[f(X)]=1N∑i=1Nf(xi)\mathbb E_{\mathcal X} [f(\mathcal X)] = \frac{1}{N} \sum_{i=1}^N f(x_i)EX[f(X)]=N1i=1∑Nf(xi)

细致平衡原则与接收率

在马尔可夫链与平稳分布中介绍过，由于马尔可夫链自身的 齐次马尔可夫假设 的性质，使得状态转移矩阵A\mathcal AA必然是一个双随机矩阵；
由于双随机矩阵的自身性质，使得只要状态转移过程足够多，各状态随机变量的概率分布必然逼近于平稳分布：
mmm表示执行了足够多次状态转移过程后达到平稳分布的状态结果;Xend\mathcal X_{end}Xend表示转移过程的终结状态。
π(Xm+1)=π(Xm+2)=⋯=π(Xend)\pi(\mathcal X_{m+1}) = \pi(\mathcal X_{m+2}) = \cdots = \pi(\mathcal X_{end})π(Xm+1)=π(Xm+2)=⋯=π(Xend)

但是，我们并不清楚：到底状态转移至什么状态，就可以得到平稳分布了？
或者可理解为：虽然平稳分布在马尔可夫链中必然存在，但是要一直等下去，等到转移次数足够多了，就大致认为其是平稳分布。

这种做法明显是不精确的。我们需要有一种方式，来判断当前时刻随机变量的概率分布是否为平稳分布？
细致平衡原则很好地满足了这样的条件，准确来说：细致平衡原则是判断当前状态随机变量Xt\mathcal X_tXt是平稳分布的充分不必要条件。
π(Xt=xi)⋅P(Xt+1=xj∣Xt=xi)=?π(Xt=xj)⋅P(Xt+1=xi∣Xt=xj)\pi(\mathcal X_{t} =x_i) \cdot P(\mathcal X_{t+1} =x_j \mid \mathcal X_t = x_i) \overset{\text{?}}{=} \pi(\mathcal X_{t} = x_j) \cdot P(\mathcal X_{t+1} =x_i \mid \mathcal X_t =x_j)π(Xt=xi)⋅P(Xt+1=xj∣Xt=xi)=?π(Xt=xj)⋅P(Xt+1=xi∣Xt=xj)
如果式子两端相等，就可以判定ttt时刻随机变量的概率分布π(Xt)\pi(\mathcal X_t)π(Xt)是平稳分布，反之同理。

即便我们已经得到了平稳分布的判定条件，但是平稳分布只是一个无限逼近的过程，仅凭判定条件，我们可能依然很难得到那个达到平稳分布的状态mmm。

而MH采样算法构建一种方式，通过细致平衡原则将逼近过程映射成概率形式。并基于该概率拒绝或者接受样本。从而在马尔可夫链的任意状态下，均可采集样本。

MH采样算法的接受率α(xi,xj)(i,j∈{1,2,⋯,K})\alpha(x_i,x_j)(i,j \in \{1,2,\cdots,\mathcal K\})α(xi,xj)(i,j∈{1,2,⋯,K})表示如下：
α(xi,xj)=min⁡[1,π(xj)P(X=xi∣X=xj)π(xi)P(X=xj∣X=xi)]\alpha(x_i,x_j) = \min \left[1, \frac{\pi(x_j)P(\mathcal X = x_i \mid \mathcal X = x_j)}{\pi(x_i)P(\mathcal X = x_j \mid \mathcal X = x_i)}\right]α(xi,xj)=min[1,π(xi)P(X=xj∣X=xi)π(xj)P(X=xi∣X=xj)]
在细致平衡中添加接收率，可以使当前状态与下一时刻状态强行细致平衡，从而将细致平衡的可能性转化为 α(xi,xj)\alpha(x_i,x_j)α(xi,xj)的接受概率：
π(Xt=xi)⋅P(Xt+1=xj∣Xt=xi)⋅α(xi,xj)=π(Xt=xj)⋅P(Xt+1=xi∣Xt=xj)⋅α(xj,xi)\pi(\mathcal X_{t} =x_i) \cdot P(\mathcal X_{t+1} =x_j \mid \mathcal X_t = x_i)\cdot \alpha(x_i,x_j) = \pi(\mathcal X_{t} = x_j) \cdot P(\mathcal X_{t+1} =x_i \mid \mathcal X_t =x_j) \cdot \alpha(x_j,x_i)π(Xt=xi)⋅P(Xt+1=xj∣Xt=xi)⋅α(xi,xj)=π(Xt=xj)⋅P(Xt+1=xi∣Xt=xj)⋅α(xj,xi)
就以α(xi,xj)\alpha(x_i,x_j)α(xi,xj)的概率采集xjx_jxj样本，在α(xi,xj)\alpha(x_i,x_j)α(xi,xj)概率结果之外，自然只能拒绝该样本，并将xix_ixi再一次加入样本集合中。

MH采样算法的弊端

基于上述流程，我们发现：

添加接受率这种操作确实能够使当前状态与下一时刻状态之间强行细致平衡，但是这种细致平衡同样是基于概率的。

如果状态转移矩阵A\mathcal AA的选择不当，导致每次状态转移后指向的样本对应的接收率都不高，此时，这个样本即便指向了它，但也大概率会被拒绝。
即出现‘下一状态对于各个结果的状态转移概率都很平均，都不太高，以至于没有明显区分’;
随机变量X\mathcal XX可能是一个高维度的随机变量，那么π(X)\pi(\mathcal X)π(X)可视作随机变量X\mathcal XX各维度分量的联合概率分布。即：
X=(x1,x2,⋯,xp)Tπ(X)=π(x1,x2,⋯,xp)\begin{aligned}\mathcal X = (x_1,x_2,\cdots,x_p)^{T} \\ \pi(\mathcal X) = \pi(x_1,x_2,\cdots,x_p) \end{aligned}X=(x1,x2,⋯,xp)Tπ(X)=π(x1,x2,⋯,xp)
实际上，基于高维样本的马尔可夫链的采样过程很困难，加上第一条问题，采出来的样本不仅量少，而且样本的有效性也很差。
这里的‘有效性差’指采出来的样本‘无法近似表示’平稳分布。

吉布斯采样方法

吉布斯采样的采样过程

吉布斯采样方法的核心思想：坐标上升法。

相比于MH算法将整个样本(所有维度)一次性采集完成，吉布斯采样使用针对某一维度进行采样，并将除采样维度外的其他维度信息进行固定；下一状态采样其他维度过程中，将上一状态采样的维度加入本次采样中进行更新，直至所有维度均采样一遍为止。

以一个基于二维信息的联合概率分布π(x1,x2)\pi(x_1,x_2)π(x1,x2)为例：

已知一个初始化样本点S=(x1(0),x2(0))\mathcal S = (x_1^{(0)},x_2^{(0)})S=(x1(0),x2(0))；
一个初始化集合D={S}\mathcal D = \{\mathcal S\}D={S}；
固定x1x_1x1维度：
- 固定x1x_1x1维度的样本信息x1(0)x_1^{(0)}x1(0)，从 后验概率P(x2∣x1(0))P(x_2 \mid x_1^{(0)})P(x2∣x1(0))中采样x2x_2x2维度的维度样本：
  x2(1)∼P(x2∣x1(0))x_2^{(1)} \sim P(x_2 \mid x_1^{(0)})x2(1)∼P(x2∣x1(0))
- 此时x2(1)x_2^{(1)}x2(1)满足 细致平衡原则：
  此时，我们可以将x2(0),x2(1)x_2^{(0)},x_2^{(1)}x2(0),x2(1)均视为‘后验概率’P(x2∣x1(0))P(x_2 \mid x_1^{(0)})P(x2∣x1(0))中随机产生的样本点~
  x2(0)→x2(1):π(x1(0),x2(0))⋅P(x2(1)∣x1(0))x2(1)→x2(0):π(x1(0),x2(1))⋅P(x2(0)∣x1(0))π(x1(0),x2(0))⋅P(x2(1)∣x1(0))=π(x1(0),x2(1))⋅P(x2(0)∣x1(0))x_2^{(0)} \to x_2^{(1)} : \quad\pi(x_1^{(0)},x_2^{(0)}) \cdot P(x_2^{(1)} \mid x_1^{(0)}) \\ x_2^{(1)} \to x_2^{(0)} : \quad \pi(x_1^{(0)},x_2^{(1)}) \cdot P(x_2^{(0)} \mid x_1^{(0)}) \\ \pi(x_1^{(0)},x_2^{(0)}) \cdot P(x_2^{(1)} \mid x_1^{(0)}) = \pi(x_1^{(0)},x_2^{(1)}) \cdot P(x_2^{(0)} \mid x_1^{(0)})x2(0)→x2(1):π(x1(0),x2(0))⋅P(x2(1)∣x1(0))x2(1)→x2(0):π(x1(0),x2(1))⋅P(x2(0)∣x1(0))π(x1(0),x2(0))⋅P(x2(1)∣x1(0))=π(x1(0),x2(1))⋅P(x2(0)∣x1(0))
- 至此，得到一个中间的过渡样本：M^=(x1(0),x2(1))\hat {\mathcal M} = (x_1^{(0)},x_2^{(1)})M^=(x1(0),x2(1))。
固定x2x_2x2维度：
- 借助样本点M^\hat {\mathcal M}M^，固定x2x_2x2维度的样本信息x2(1)x_2^{(1)}x2(1)，从 后验概率P(x1∣x2(1))P(x_1 \mid x_2^{(1)})P(x1∣x2(1))中采样x1x_1x1维度的维度样本：
  x1(1)∼P(x1∣x2(1))x_1^{(1)} \sim P(x_1 \mid x_2^{(1)})x1(1)∼P(x1∣x2(1))
- 此时x1(1)x_1^{(1)}x1(1)满足 细致平衡原则：
  同理，此时我们可以将x1(0),x1(1)x_1^{(0)},x_1^{(1)}x1(0),x1(1)均视为‘后验概率’P(x1∣x2(1))P(x_1 \mid x_2^{(1)})P(x1∣x2(1))中随机产生的样本点~
  π(x1(0),x2(1))⋅P(x1(1)∣x2(1))=π(x1(1),x2(1))⋅P(x1(0)∣x2(1))\pi(x_1^{(0)},x_2^{(1)}) \cdot P(x_1^{(1)} \mid x_2^{(1)}) = \pi(x_1^{(1)},x_2^{(1)}) \cdot P(x_1^{(0)} \mid x_2^{(1)})π(x1(0),x2(1))⋅P(x1(1)∣x2(1))=π(x1(1),x2(1))⋅P(x1(0)∣x2(1))
- 至此，再次得到一个样本：M=(x1(1),x2(1))\mathcal M = (x_1^{(1)},x_2^{(1)})M=(x1(1),x2(1))。
注意，至此我们将所有维度全部固定一次(x1,x2x_1,x_2x1,x2)。并得到一个新的样本点M=(x1(1),x2(1))\mathcal M = (x_1^{(1)},x_2^{(1)})M=(x1(1),x2(1))。将该样本点添加到样本集合D\mathcal DD中：
D={S,M}\mathcal D = \{\mathcal S,\mathcal M\}D={S,M}
后续步骤将以M\mathcal MM点为基准点，重复执行上述操作，采样出新的样本点，直至样本集合D\mathcal DD足够大，最终进行蒙特卡洛方法求解。
图形描述表示如下：

吉布斯采样的推导过程

为什么可以使用这种迭代方式进行采样？这种采样方式本质上是MH算法中的一种特殊情况：接受率为1的MH采样方法。

为什么这种采样方式得到的样本都可以完整接受它？

构建如下场景：
假设要采样的隐变量Z\mathcal ZZ是一个ppp维的随机变量：
Z=(z1,z2,⋯,zp)T\mathcal Z = (z_1,z_2,\cdots,z_p)^{T}Z=(z1,z2,⋯,zp)T
该随机变量Z\mathcal ZZ对应在各时间状态中的概率分布π(Z)\pi(\mathcal Z)π(Z)可看作是关于Z\mathcal ZZ的联合概率分布：
π(Z)=π(z1,z2,⋯,zp)\pi(\mathcal Z) = \pi(z_1,z_2,\cdots,z_p)π(Z)=π(z1,z2,⋯,zp)
针对上述吉布斯采样的采样过程，任意维度的样本zi∗(i=1,2,⋯,p)z_i^*(i=1,2,\cdots,p)zi∗(i=1,2,⋯,p)可视作从状态转移矩阵P(zi∣z1,⋯,zi−1,zi+1,⋯,zp)P(z_i \mid z_1,\cdots,z_{i-1},z_{i+1},\cdots,z_p)P(zi∣z1,⋯,zi−1,zi+1,⋯,zp)中生成的一个样本：
zi∗∼P(zi∣z1,⋯,zi−1,zi+1,⋯,zp)z_i^* \sim P(z_i \mid z_1,\cdots,z_{i-1},z_{i+1},\cdots,z_p)zi∗∼P(zi∣z1,⋯,zi−1,zi+1,⋯,zp)
同理，在生成其他维度样本过程中，如zj∗(j≠i)z_j^*(j \neq i)zj∗(j=i)，需要将前面步骤采集的zi∗z_i^*zi∗样本加入到后验概率分布中：
zj∗∼P(zj∣z1,⋯,zi∗,⋯,zj−1,zj+1,⋯,zp)z_j^* \sim P(z_j \mid z_1,\cdots,z_i^*,\cdots,z_{j-1},z_{j+1},\cdots,z_p)zj∗∼P(zj∣z1,⋯,zi∗,⋯,zj−1,zj+1,⋯,zp)
其他维度样本均同理。直至所有维度样本均采样一次后，得到的结果是一个新的样本点。

上述过程中，吉布斯采样存在一个特殊现象：对于一个完整样本的ppp次迭代中，每一次迭代转移概率都是固定的。因此即便存在维度替换的情况，但并不影响当前步骤的后验概率结果。

示例：
从如下两个样本中采样维度jjj的样本zjz_jzj：
zj∼P(zj∣z1,⋯,zi∗,⋯,zj−1,zj+1,⋯,zp)zj∼P(zj∣z1,⋯,zi,⋯,zj−1,zj+1,⋯,zp)z_j \sim P(z_j \mid z_1,\cdots,z_i^*,\cdots,z_{j-1},z_{j+1},\cdots,z_p) \\ z_j \sim P(z_j \mid z_1,\cdots,z_i,\cdots,z_{j-1},z_{j+1},\cdots,z_p)zj∼P(zj∣z1,⋯,zi∗,⋯,zj−1,zj+1,⋯,zp)zj∼P(zj∣z1,⋯,zi,⋯,zj−1,zj+1,⋯,zp)

这两个分布的唯一区别在于：第iii维度是否被替换过。但不需要担心这种情况。因为上述两次采样绝不会出现在同一个样本的迭代过程中。

观察接受率中的分数项，并表示为吉布斯采样的概率形式：
注意，此时的Z\mathcal ZZ已经是一个高维向量;z−iz_{-i}z−i表示除去ziz_izi分量后的剩余分量。
P(Z∗∣X)⋅Q(Z∣Z∗)P(Z∣X)⋅Q(Z∗∣Z)=π(zi∗∣z−i∗)⋅π(z−i∗)⋅p(zi∣z−i∗)π(zi∣z−i)⋅π(z−i∗)⋅p(zi∗∣z−i)=1\begin{aligned} & \frac{P(\mathcal Z^* \mid \mathcal X) \cdot \mathcal Q(\mathcal Z \mid \mathcal Z^*)}{P(\mathcal Z \mid \mathcal X) \cdot \mathcal Q(\mathcal Z ^* \mid \mathcal Z)} \\ & = \frac{\pi(z_i^* \mid z_{-i}^*) \cdot \pi(z_{-i}^*)\cdot p(z_i \mid z_{-i}^*)}{\pi(z_i \mid z_{-i}) \cdot \pi(z_{-i}^*) \cdot p(z_i^* \mid z_{-i})} \\ & = 1 \end{aligned}P(Z∣X)⋅Q(Z∗∣Z)P(Z∗∣X)⋅Q(Z∣Z∗)=π(zi∣z−i)⋅π(z−i∗)⋅p(zi∗∣z−i)π(zi∗∣z−i∗)⋅π(z−i∗)⋅p(zi∣z−i∗)=1

吉布斯采样的代码实现

基于上述的吉布斯采样过程，我们通过采样的方式近似一个给定的2维高斯分布：

首先，这个近似操作只是一个简化了的工作：

二维高斯分布它的各维度也是高斯分布；
各维度之间相互独立，因此，状态转移概率(条件概率)就是边缘概率分布本身；

说白了，就是从两个独立的高斯分布中进行采样。但是吉布斯采样的过程还是要保留的。
代码表示如下：

import math
import matplotlib.pyplot as plt
import randomdef one_dim_pdf(x,mu,sigma):return (1 / (math.sqrt(2 * math.pi) * sigma)) * math.exp(-1 * (((x - mu) ** 2)/ (2 * (sigma ** 2))))def get_sample(mu,sigma):while True:random_value = random.uniform(-5, 15)u = random.uniform(0, 1)prob = one_dim_pdf(random_value, mu, sigma)if u < prob:breakreturn random_valuedef gibbs(start_point,mu_1,mu_2,sigma_1,sigma_2,loop_num):assert len(start_point) == 2start = start_pointfor i in range(loop_num):start_x, start_y = start[0], start[1]randomv = get_sample(mu_1,sigma_1)plt.scatter([start_x,start_x],[start_y,randomv],c="k",s=2)randomq = get_sample(mu_2,sigma_2)plt.scatter([start_x,randomq],[randomv,randomv],c="k",s=2)start = [randomq,randomv]plt.show()if __name__ == '__main__':mu_1 = 0sigma_1 = 1mu_2 = 10sigma_2 = 2start_point = [0,1]gibbs(start_point,mu_1,mu_2,sigma_1,sigma_2,loop_num=500)

对应图像表示如下：

下一节将介绍马尔可夫链蒙特卡洛方法面对的挑战

相关参考：
机器学习-白板推导系列-蒙特卡洛方法4（Monte Carlo Method）-吉布斯采样（Gibbs）

机器学习笔记之马尔可夫链蒙特卡洛方法(四)吉布斯采样相关推荐

机器学习笔记之马尔可夫链蒙特卡洛方法(三)MH采样算法
机器学习笔记之马尔可夫链蒙特卡洛方法--MH采样算法引言回顾:马尔可夫链与平稳分布马尔可夫链平稳分布 MH采样算法采样思路 MH采样算法过程引言上一节介绍了马尔可夫链(Markov Ch ...
机器学习笔记马尔可夫链蒙特卡洛方法(二)马尔可夫链与平稳分布
机器学习笔记之马尔可夫链蒙特卡洛方法--马尔可夫链与平稳分布引言回顾:蒙特卡洛方法马尔可夫链与平稳分布马尔可夫链平稳分布细致平衡关于平稳分布的补充马尔可夫链的本质平稳分布的收敛性证明 ...
不用数学也能讲清贝叶斯理论的马尔可夫链蒙特卡洛方法？这篇文章做到了
作者 | Ben Shaver 翻译 | 刘畅编辑 | Donna 大多数时候,贝叶斯统计在结果在最好的情况下是魔法,在最糟糕时是一种完全主观的废话.在用到贝叶斯方法的理论体系中,马尔可夫链蒙特卡洛 ...
不用数学讲清马尔可夫链蒙特卡洛方法？
翻译 | 刘畅大多数时候,贝叶斯统计在结果在最好的情况下是魔法,在最糟糕时是一种完全主观的废话.在用到贝叶斯方法的理论体系中,马尔可夫链蒙特卡洛方法尤其神秘. 这篇文章将介绍马尔可夫链蒙特卡洛方法, ...
AI算法连载04：数学基础之蒙特卡洛方法与MCMC采样
相关文章: AI算法连载01:数学基础之线性代数 AI算法连载02:概率论基础 AI算法连载03:数学基础之数值计算导语:不懂算法的工程师做不了AI,不懂算法的产品经理将把AI带入泥潭. 一.前言 ...
一份数学小白也能读懂的「马尔可夫链蒙特卡洛方法」入门指南
文章来源:机器之心在众多经典的贝叶斯方法中,马尔可夫链蒙特卡洛(MCMC)由于包含大量数学知识,且计算量很大,而显得格外特别.本文反其道而行之,试图通过通俗易懂且不包含数学语言的方法,帮助读者对 M ...
李宏毅机器学习笔记(十六)——无监督学习(四)：自编码器
文章目录一.自编码器的引入二.自编码器的应用 1.文本检索 2.寻找相似图片 3.预训练神经网络 4.面向图片的自编码器 5.进行内容生成三.自编码器的改进 1.去除噪音 2.重定义损失函数 3 ...
什么是吉布斯采样（Gibbs Sampling）
目录 1 蒙特卡洛方法 1.1 蒙特卡洛方法的作用 1.2 非均匀分布采样 1.3 分布p(x)不好采样怎么办? 2 什么是吉布斯采样 2.1 马尔可夫链 2.1.1 什么是马尔可夫链呢? 2.1.2 ...
【机器学习基础】数学推导+纯Python实现机器学习算法21：马尔可夫链蒙特卡洛...
Python机器学习算法实现 Author:louwill Machine Learning Lab 蒙特卡洛(Monte Carlo,MC)方法作为一种统计模拟和近似计算方法,是一种通过对概率模型随 ...

机器学习笔记之马尔可夫链蒙特卡洛方法(四)吉布斯采样

机器学习笔记之马尔可夫链蒙特卡洛方法——吉布斯采样

引言