采样方法 | 记录思考

在前两周组内的技术分享中，分享了采样方法。之前在研究生阶段就对采样方法很是疑惑，特别是在看LDA时，用到的Gibbs采样，很多次尝试去学习这一知识点，但都一知半解。所以，借这个机会认真学习一下采样方法相关的知识。本文主要是记录一下自己在学习采样方法时，对不同采样方法原理的理解，主要包括蒙特卡洛方法介绍和5中不同的采样方法。

蒙特卡洛方法

首先说一下蒙特卡洛方法，Monte Carlo方法，又称为统计模拟法、随机抽样技术，是一种随机模型方法，以概率和统计理论为基础的一种计算方法。是使用随机数（或伪随机数）来解决很多复杂问题的计算方法。其核心就是通过将所要求解的问题同一定的概率模型相联系，利用计算机进行模拟或抽样，以获得问题的近似解。简单地理解蒙特卡洛方法，其实就是将要解决的采用一定的方式进行转换，转换之后的问题通过利用计算机实现统计模拟或抽样，从而获得问题是近似解。至于为什么叫Monte Carlo方法，可能是和闻名世界的赌城——摩纳哥的一个小山城有关吧！具体是什么原因，不必去深究。

来看一个典型的用Monte Carlo方法解决实际问题的例子。计算圆周率Π的值，如果采用Monte Carlo方法进行计算，可以进行如下的转换：

首先，设置一个边长为1的正方形，其内有一圆，圆的半径为0.5，如下图所示：
随机的向正方形内打点，打的点在圆内的概率等于圆与正方形的面积之比0.25Π
随机产生M个点(x,y)，其中，x和y都是区间[0,1]内的符合均匀分布的随机数
假设落在圆内的点有N个，当M足够大时，根据大数定理，频率等于概率，就有：$$\frac {N} {M} = 0.25Π$$，$$Π=\frac {4N} {M}$$

实际用代码去模拟一下，得出的圆周率的值为：3.143748，和3.1415926非常接近。

def samplePI(maxCnt = 1000000):
accCnt = 0
for i in range(maxCnt):
	x = np.random.random()
	y = np.random.random()
	if np.sqrt(x**2 + y**2) < 1:
		accCnt += 1
print "PI=", float(accCnt) / maxCnt *4

通过一个简单的例子，可以直观地理解Monte Carlo方法。采样Monte Carlo方法去解决实际问题一般可以分为如下的步骤：

对复杂的问题进行转换，构造或描述随机过程。例如，将计算圆周率的问题，转换为在一个正方形内进行打点，将圆周率的值和点在圆内的概率联系起来。
从已知的概率分布中进行采样,采样出符合特定分布的样本。例如，从符合[0,1]均匀分布中采样出x和y。
建立估计量进行计算。

Monte Carlo方法是一种通用的计算技术，可以解决如下问题:

随机模拟：从一个pdf产生”典型”的样本
计算积分：在高维空间中的积分
优化问题
学习：MLE:f(x;Θ)

利用Monte Carlo方法去解决实际问题的第2步是从符合特定分布中采样出样本。我们知道，计算机本身无法产生真正的随机数，但可以根据一定的算法产生伪随机数。比如，通过线性同余发生器可以生成伪随机数，我们可以用确定性的算法生成[0,1]之间符合均匀分布的随机数，可以被当成真正的随机数使用。而至于产生符合比较复杂分布的随机数，则是以均匀分布为基础进行采样，以获得符合特定复杂分布的样本，这就是采样方法。不同的采样方法，采样过程不同，但大都基于均匀分布来进行的。下面就分别介绍，几种不同的采样方法。

Inverse Transform Sampling

在介绍不同的采样方法之前，先说一下概率分布，概率分布一般分为连续型和离散型两种，离散型的概率分布一般用概率质量分布函数(pmf)表示，而连续分布用概率密度(pdf)表示,对pmf进行累加或对pdf进行积分的函数，对应于累积分布函数(cdf)。所有pmf的取值之和为1，对pdf在其定义域上进行积分，积分的值为1。
对于一些简单的分布p(x),我们可以直接进行采样，比如，采样符合p(x)=[0.1,0.3,0.5,0.1]分布的样本 $x=x_1,x_2,x_3,x_4$分布，我们可以直接从Uniform(0,1)采样出一个数，x小于0.1，则采样的值为$x_1$,x等于0.1小于0.4则为$x_2$，以此进行类推。对于比较复杂的分布p(x)，不能直接进行采样。如果要采样的目标分布为p(x),它的累积分布函数F(x)能够求出来，F(x)的反函数也能求出来，那么，Inverse Transform Sampling的采样过程如下：

1
2
3

Inverse Transform Sampling
1. 从Uniform(0,1)中随机采样一个点，用u表示。
2. 计算$F^-1(u)$的值为x,则x就是服从p(x)的分布的一个样本点。

简单的证明一下：
设F(x)为目标采样分布P(x)的累积分布，$x=F^{-1}(u) u\in[0,1]$为F(x)的反函数。因为F(x)是单调递增的(累积函数的性质)，所以$x=F^{-1}(u)$也是单调递增函数，对于如下不等式:
$$F^{-1}(u)\leq F^{-1}(F(x)), if u \leq F(x) （1）$$
根据反函数的定义有：
$$F^{-1}(u) \leq x, if u \leq F(x) （2）$$
根据Uniform(0,1)的定义，其累积分布函数如下：

所以，采样出的点的$F^{-1}(u)$累积分布函数如下：
$$P(F^{-1}(u) \leq x)=P(u \leq F(x))=H(F(x))=F(x)$$由此可见，采样出的点的累积分布函数为F(x),所以，为符合p(x)分布的样本。
虽然，利用这种方法可以采样出符合特定分布的样本，但是这种采样方法存在一个缺陷，即，如果要采样的分布p(x)累积分布函数不能够求出来或者累积分布函数没有反函数，这种方法就失效了

Acceptance-Rejection Sampling

对于采样的目标分布p(x),如果其比较复杂，可以采用Acceptance-Rejection Sampling接受-拒绝采样，来进行采样。这种采样方法不直接对p(x)进行采样，而是选择另一个分布q(x)，存在常数M，使得$p(x)\leq Mq(x),\forall x $,称为建议分布(Proposal Density),这个分布容易进行采样，通过对q(x)的采样，来实现对p(x)的采样。其具体的采样过程如下：

Acceptance-Rejection Sampling的采样过程：
1. 从q(x)中随机采样出一个样本点x0
2. 在从Uniform(0, Mq(x0))中采样出另一个样本点u
3. 如果u>p(x0),则拒绝x0,并且重复前面的步骤。否则接受x0为符合p(x)分布的样本

这种采样方法可以用计算圆周率的方法来进行理解。把正方形看做一个分布Mq(x),圆形看做要采样的目标分布q(x),通过对Mq(x)进行随机采样，即打点，如点果落在圆内，表明该点符合圆形这个分布，即可以做符合p(x)的样本。即，如下图所示的圆内的点都是符合圆形这个分布的点。

Acceptance-Rejection Sampling，通过对q(x)的采样，实现了对p(x)的采样。这种采样方法的接受率正比于$\frac {1} {M}$,等于p(x)下面的面积除以Mq(x)下面的面积。可以看出，只有当M尽可能的小时，采样出的点被接受的概率才会大，从而可以提升采样的效率，因此在使用该方法进行采样是M的选择比较重要。同时,可以看出，如果q(x)和p(x)的形相似时，M就会越小，接受率就会越高。然而，对于高维的目标采样分布，q(x)可能不容易寻找，且M也可能会很大，此时，接受率就会变小，采样效率会变差。

Importance Sampling

Importance Sampling 这种采样方法，其并不是为了获得符合特定分布p(x)的样本，而是为了解决当p(x)不容易进行采样时，计算E[f(x)]，x符合p(x)分布，的问题。其可以进行如下的转换：
首先，根据期望的定义，E[f(x)]的计算公式如下：
$$
E[f(x)] = \int_x f(x)p(x)dx
$$
因为，p(x)的样本不容易获取，Importance Sampling同样引入一个建议分布q(x)，比较容易获取符合q(x)分布的样本，进行如下转换（以下内容参考：http://blog.csdn.net/dark_scope/article/details/70992266）：
$$
\int_x f(x)p(x)dx = \int_x f(x)\frac {p(x)} {q(x)} q(x) dx
= \int_x g(x)q(x)dx where g(x)=f(x) \frac {p(x)} {q(x)} = f(x)w(x)
$$
可以看出，通过上面的转化，可以将计算f(x)在p(x)下的期望转化为求g(x)在q(x)分布下的期望。其中，$w(x) = \frac {p(x)} {q(x)}$被称为Importance Weight。
但是，有些时候p(x)也是很难计算的，更常见的情况是比较方便的计算$\hat p (x)$和$\hat q(x)$
$$p(x) = \frac {\hat p(x) }{Z_p}$$$$p(x) = \frac {\hat p(x)}{Z_p}$$
其中，$Z_{p/q}$是一个标准化项，可以看成是一个常数，是的$\hat p(x)$或者$\hat q(x)$等比例变化为一个概率分布，也可以理解为softmax里的分母。其中：
$$Z_p=\int_x \hat p(x)dx$$$$Z_q=\int_x \hat q(x)dx$$
在这种情况下，Importance Sampling 可以进行如下的转换：
$$ \int_x f(x)p(x)dx=\int_x f(x)\frac {p(x)} {q(x)} q(x)dx\\\\ =\int_x f(x) \frac {\hat p(x)/Z_p} {\hat q(x)/Z_q} q(x)dx\\\\=\frac {Z_q}{Z_p} \int_x f(x) \frac {\hat p(x)}{\hat q(x)}\\\\=\frac {Z_q}{Z_p} \int_x \hat g(x)q(x)dx\\\\其中，\hat g(x) = f(x)\frac {\hat g(x)}{\hat q(x)}=f(x)\hat w(x)
$$
而$\frac {Z_q}{Z_p}$直接计算不太好计算，而它的倒数：
$$\frac {Z_p}{Z_q}=\frac {1}{Z_q}\int_x \hat p(x)dx, Z_q=\frac {\hat q(x)} {q(x)} $$
所以：
$$
\frac {Z_p}{Z_q}=\frac {\hat p(x)}{\hat q(x)}q(x)dx = \int_x \hat w(x)q(x)dx
$$
这样，假设能方便从q(x)进行采样，所以上式又被转换为一个Monte Carlo可解的问题，也就是说：
$$
\frac {Z_p}{Z_q}=\frac {1}{m} \sum_{i=1}^m \hat w(x_i), x_i符合q(x)分布。
$$
最终，求解E[f(x)]的问题可以转换为：
$$E[f(x)]=\frac {1}{m} \sum_{i=1}^m \hat w(x_i)f(x_i), 其中，x_i为符合q(x)的样本\\\ \hat w(x_i)=\frac {\hat w(x_i)}{\sum_{i=1}^m \hat w(x_i)}$$
所以，我们可以在不用知道q(x)确切值的情况下，就可以近似地计算得到E[f(x)]。其计算过程如下：
Importance Sampling采样过程：

首先为p(x)找到一个建议分布q(x),q(x)比较容易采样。
然后从q(x)中采样出m个点x
带入$E[f(x)] = \frac {1}{m} \sum_{i=1}^m \hat w(x_i)f(x_i)$ 计算期望。
虽然这种方法能够work,但是在高维空间里找到一个这样合适的q(x)非常难。即使有 Adaptive importance sampling 和 Sampling-Importance-Resampling的出现，要找到一个同时满足容易抽样并且和目标分布相似的建议分布，通常是不可能的！

MCMC: Markov Chain Monte Carlo

Importance Sampling和Acceptance-Rejection Sampling虽然能够实现对一些分布的采样，但是只有当选取的建议分布q(x)和要进行采样的目标分布p(x)很近似时才表现好，所以选取合适的q(x)是非常关键。当在高维空间进行采样，标准的采样方法会失败，对于Acceptance-Rejection Sampling,当目标分布的维数增高时，拒绝率会趋近于100%，采样的效率会很低。对于Importance Sampling，大多数的样本的权重值会趋近于0。对于高维复杂问题，可以用马尔科夫链（Markov Chain)产生一系列相关样本，实现对目标分布的采样。MCMC是一种用一定范围内的均匀分布的随机数对高维空间概率进行采样的通用技术，其基本思想是设计一个马尔科夫链，使得其稳定概率分布为要采样的目标分布$\pi(x)$

首先来看一下马尔科夫链的定义及其平稳分布：

马尔科夫性质：某一时刻状态转移的概率只依赖于它前一个状态
定义：假设存在状态序列$… X_{t-2},X_{t-1},X_t,X_{t+1}…$,时刻t+1的状态的条件概率只依赖于t时刻的状态$x_t$,即：
$$P(X_{t+1}|…X_{t-2},X_{t-1},X_t)=P(X_{t+1}|X_t)$$
马尔科夫链：满足马尔科夫性质的随机过程

以天气变化来解释一下上面的定义，假设每天的天气是一个状态的话，状态转移可以看成是天气的变化，比如从晴天变成阴天、从阴天变成雨天等。马尔科夫性质讲的是，今天的天气情况只依赖于昨天的天气，和前天以及之前的天气状况没有任何关系。马尔科夫链可以看成每天的天气按照这个规律进行变化的一个过程。一个马尔科夫链可以由下面的公式定义：

一个马尔科夫链一般由三部分构成：

状态空间：可以理解为天气状况的所有情况 {阴天，晴天，雨天，…}
初始状态：可以理解为第一天的天气情况
状态转移矩阵：可以理解为所有由一种天气状况变为另一种天气状况的概率

马尔科夫链具有一个非常重要的性质：马尔科夫链的平稳分布。来看一个例子，假设一个国家的人口地域分布分为：农村、城镇和城市3种状态。每年人口流动情况如下图：

上图表示每种状态转移到另一种状态的概率。如果定义矩阵P，P的某一位置P(i,j)的值为P(j|i)，表示从状态i转换为状态j的概率，则根据上图可以得到马尔科夫链的状态转移矩阵为：

假设初始状态的人口地域分布为$\pi_0=[\pi_0(1),\pi_0(2),\pi_0(3)]$，每年人口按照状态转移矩阵P进行转移，n年后人口的地域分布为$\pi_n=\pi_{n-1}P=…=\pi P^n$。假设存在如下两种初始的人口分布:

$\pi_0=[0.5,0.4,0.1]$
$\pi_0=[0.3,0.4,0.3]$

按照状态转移矩阵进行转移一定年数后的人口分布情况分布如下图所示：

可以看出，尽管采用了不同的初始化状态，但最终的概率分布都趋近于一个稳定的概率分布[0.167,0.388,0.444]。可以看到，马尔科夫链模型的状态转移矩阵收敛到稳定概率分布和初始状态概率分布无关。也就是说，如果得到了稳定概率分布对应的马尔科夫链模型的状态转移矩阵，我们可以从任意的概率分布样本开始，带入马尔科夫链模型的状态转移矩阵，经过一系列的状态转移，最终样本的分布会趋近于稳定的概率分布。用数学的语言来定义一下马尔科夫链的收敛性质：
如果一个非周期的马尔科夫链，其状态转移矩阵为P，并且它的任何两个状态之间是联通的，那么$\lim_{n \rightarrow +\infty} P_{ij}^n$ 与i无关，则有：

$\lim_{n \rightarrow +\infty} P^n = \begin{bmatrix}
{\pi(1)}&{\pi(2)}&{\cdots}&{\pi(n)}\\\\
{\pi(1)}&{\pi(2)}&{\cdots}&{\pi(n)}\\\\
{\vdots}&{\vdots}&{\ddots}&{\vdots}\\\\
{\pi(1)}&{\pi(2)}&{\cdots}&{\pi(n)}\\\\
\end{bmatrix}$
$\pi(j)=\sum_{i=0}^\infty P_{ij}$，$p_{ij}$表示状态i转移到状态j的概率
$\pi$是方程$\pi P = \pi$的唯一非负解，其中，$$\pi=[\pi(1),\pi(2),….,\pi(j),…] \sum_{i=0}^\infty \pi(i)=1$$

上面的性质中有如下的说明：

非周期的马尔科夫链：主要是指马尔科夫链的状态转化不是循环的，如果是循环的则永远不会收敛。我们一般遇到的马尔科夫链一般都是非周期性的。用数学方式可以表述为：对于任意某一状态i,d为集合${n | n\geq 1,P_{ij}^n > 0}$的最大公约数为1，如果d=1，则该状态为非周期的。
任何两个状态是联通的：指从任意一个状态可以通过有限步到达其他的任意一个状态，不会出现条件概率为0导致不可达的情况。
马尔科夫链的状态可以是有限的，也可以是无限的。因此，可以用于连续概率分布和离散概率分布。
$\pi$通常称为马尔科夫链的平稳分布。

基于马尔科夫链的采样方法

从马尔科夫链的收敛性质可以看出，如果我们能够得到某个平稳分布对应的马尔科夫链状态转移矩阵，就很容易采样出这个平稳分布的样本。
假设任意初始化的概率分布为$\pi_0(x)$，经过第一轮转移后的概率分布是$\pi_0(x)$，经过i轮后的概率分布为$\pi_i(x)$。假设经过n轮后马尔科夫链收敛到平稳分布$\pi(x)$,即：
$$\pi_n(x)=\pi_{n+1}(x)=…=\pi(x)$$
那么经过n轮之后的，沿着状态转移矩阵进行转移得到的样本都是符合$\pi(x)$分布的样本。也就是说，如果从一个具体的初始状态$x_0$出发，沿着马尔科夫链按照状态转移矩阵进行跳转，假设n次跳转后收敛，那么得到一个转移状态序列$X_0,X-1,…,X_n,X_{n+1},….$，则$X_n,X_{n+1},…$都是符合平稳分布$\pi(x)$的样本。

基于马尔科夫链的采样过程如下：

输入：状态转移矩阵P，设定状态转移次数$n_1$，需要采样的样本个数$n_2$
从任意概率分布采样得到一个初始状态值$x_0$
for t=0 to $n_1+n_2 - 1$：从条件概率分布$P(x|x_t)$中采样出样本为$x_{t+1}$
输出：样本$(x_{n_1},x_{n_1+1},…,x_{n_1+n_2})$即为符合平稳分布$\pi(x)$的样本。

解释一下上面的从条件概率分布$P(x|x_t)$中采样出样本为$x_{t+1}$。为什么要从$P(x|x_t)$中进行采样？
因为，当前的状态为$x_t$，从当前状态$x_t$进行转移，可转移到的状态为状态空间集合，所以要以当前状态为条件向其他状态进行转移，所以要从$P(x|x_t)$中进行采样。例如，人口分布的例子，假设当前的状态是农村，那么从农村转移到农村、城镇和城市的概率分别为0.5、0.4和0.1，即$P(x|x_t)=[0.5,0.4,0.1]$，转移时要按照这个分布进行转移，即采样时要按照这个分布进行采样，使用均匀分布可以很容易实现。而状态转移为连续的情况，则$P(x|x_t)$则是一个具体的连续分布，通过对$P(x|x_t)$的采样，实现转移。
可以看出，我们要采样的目标分布为$\pi(x)$,如果我们能够构建一个状态转移矩阵P，使得其马氏链上的平稳分布是$\pi(x)$，那么初始化一个状态，按照状态转移矩阵P进行转移，经过n次后收敛，第n+1次后的样本都是符合$\pi(x)$分布的。所以，MCMC采样方法的关键是构建状态转移矩阵P。如何构建这个状态转移矩阵P呢？首先看一下马尔科夫链的细致平稳条件：
如果一个非周期的马尔科夫链状态转移矩阵为P和概率分布$\pi(x)$，对于所有的i，j满足：$$
\pi(i)P(i,j)=\pi(j)P(j,i)$$其中，P(i,j)表示状态i转移到状态j的概率,则称$\pi(x)$是状态转移矩阵P的平稳分布。
证明：$$
\sum_{i=1}^\infty \pi(i)P(i,j)=\sum_{i=1}^\infty \pi(j)P(j,i)=\pi(j)\sum_{i=1}^\infty P(j,i) = \pi(j)
$$写成矩阵的形式，即：$$\pi P=\pi$$由于$\pi$是$\pi P=\pi$的解，所以$\pi$是一个平稳分布。上式被称为细致平稳条件(detailed balance condition)。
其实这个定理是显而易见的，因为细致平稳条件的物理含义就是对于任何两个状态i,j 从i转移出去到j而丢失的概率质量，恰好会被从j转移回i的概率质量补充回来，所以状态i上的概率质量$\pi(i)$是稳定的，从而$\pi(x)$是马尔科夫链的平稳分布。从细致平稳条件可以看出，只要找到可以使概率分布$\pi(x)$满足细致平稳分布的矩阵P即可。
假设在进行采样之前已经存在一个状态转移矩阵Q，Q(i,j)表示从状态i转移到状态q的概率，也可以表示为Q(j|i),通常情况下：
$$
\pi(i)Q(i,j)\neq \pi(j)Q(j,i)
$$
也就是说不满足细致平稳条件。不过可以对其进行改造，使其满足细致平稳条件。具体的是引入$\alpha$使得上面的公式成立，即：
$$\pi(i)Q(i,j)\alpha(i,j)=\pi(j)Q(j,i)\alpha(j,i)$$
什么样的$\alpha$能够使上式成立呢？最简单的，按照对称性，可以取如下的值：
$$
\alpha(i,j)=\pi(j)Q(j,i)\\\\
\alpha(j,i)=\pi(i)Q(j,i)
$$
在改造Q的过程中，引入的$\alpha(i,j)$称为接受率，取值在[0,1]之间，物理意义可以理解为在原来的马尔科夫链上，从状态i以Q(i,j)的概率转跳转到状态j的时候，我们以$\alpha(i,j)$的概率接受这个转移。这样就可以得到新的马尔科夫链的转移概率为$Q(i,j)\alpha(i,j)$。通过这种改造，就能够进行采样了。MCMC采样过程如下：
MCMC采样算法

初始化马尔科夫链初始状态$X_0=x_0$
fo t=0,1,2,… do:
2.1 在t时刻马尔科夫链状态为$X_t=x_t$，从$Q(x|x_t)$中采样出一个样本y
2.2 从均匀分布中采样出u~Uniform(0,1)
2.3 如果$u<\alpha(x_t,y)=\pi(y)Q(x_t|y)$，则接受转移$x_{t+1}=y$
2.4 否则，不接受转移，$x_{t+1}=x_t$

上面的MCMC采样算法已经能够很好的工作了，但是它存在一个缺陷:马氏链在转移的过程中的接受率$\alpha(i,j)$可能偏小，这样采样过程中，马氏链不容易转移，一直处于原地，导致收敛到平稳分布的速度偏慢。假设$\alpha(i,j)=0.1,\alpha(j,i)=0.2$，假设在此时满足细致平稳条件，于是有：
$$\pi(i)Q(i,j).0.1 = \pi(j)\alpha(j,i).0.2$$将上式两边同时扩大5倍，等式变为：
$$\pi(i)Q(i,j).0.5 = \pi(j)\alpha(j,i).1$$
可以看到细致平稳条件并没有被打破，而接受率变大了。因此，我们可以把细致平稳条件中的$\alpha(i,j)$和$\alpha(j,i)$同比例放大，使得两个数中较大的那个放大到1，这样就提高了采样中的跳转的接受率。于是$\alpha(i,j)$可以取下面的值:$$
\alpha(i,j) = min{\frac {\pi(j)Q(j,i)} {\pi(i)Q(i,j)},1}
$$
这样，就完成了对MCMC采样算法的改造，这就是Metropolis-Hastings采样算法，其采样过程如下：
Metropolis-Hastings采样算法

初始化马尔科夫链初始状态$X_0=x_0$
for t=0,1,2,… do:
2.1 在t时刻马尔科夫链状态为$X_t=x_t$，从$Q(x|x_t)$中采样出一个样本y
2.2 从均匀分布中采样出u~Uniform(0,1)
2.3 如果$u<\alpha(x_t,y)=min{\frac {\pi(y)Q(y,x_t)} {\pi(x_t)Q(x_t,y)},1}$，则接受转移$x_{t+1}=y$
2.4 否则，不接受转移，$x_{t+1}=x_t$

Gibbs Sampling

虽然MCMC采样和Metropolis-Hasting采样算法已经能够解决蒙特卡罗方法中需要的任意概率分布的样本的问题。但是还是存在一定的缺陷，首先是采样过程中要计算接受率，在高维时，计算量大，可能存在辛辛苦苦计算出的接受率，最终被拒绝，不跳转。并且由于接受率的原因导致算法的收敛时间变长。其次是，对于高维空间，状态的条件概率分布好求解，但是联合分布不好求。针对这一问题，对于高维空间的数据采样，Stuart Geman和Donald Geman这两兄弟于1984年提出来了Gibbs Sampling算法。Gibbs Sampling算法思想是通过构建状态转移矩阵，使得接受率为1，从而提升了接受的效率。对于二维情形，假设存在一个概率分布P(x,y)，对于x坐标相同的两个点$A(x_1,y_1),B(x_1,y_2)$可以得到如下公式：
$$
p(x_1,y_1)p(y_2|x_1) = p(x_1)p(y_1|x_1)p(y_2|x_1)\\\\
p(x_1,y2)p(y_1|x_1)=p(x_1)p(y_2|x_1)p(y_1|x_1)
$$
上面的转换是利用乘法公式进行转换的。于是可以得到：
$$
p(x_1,y_1)p(y_2|x_1) = p(x_1,y_2)p(y_1,x_1)\\\\
p(A)p(y_2|x_1) = p(B)p(y_1|x_1)
$$
可以看到，在$x=x_1$这条直线上，如果使用条件概率分布$p(y|x_1)$作为任意两点转移概率，那么任何两点之间的转移概率满足细致平稳条件。同样，在$y=y_1$这条直线上任取两点$A(x_1,y_1),C(x_2,y_1)$，同样有：
$$p(A)p(x_2|y_1)=p(C)p(x_1|y_1)$$
于是可以构造平面上任意两点之间的转移概率矩阵Q:
$$
Q(A->B)=p(y_B|x_1), if x_A=x_B=x_1\\\\
Q(A->C)=P(x_C|y_1), if y_A=y_C=y_1\\\\
Q(A->D)=0，其他
$$
则对于平面上任意两点X,Y,很容易验证是否满足细致平稳条件：
$$p(X)Q(X->Y)=P(Y)Q(Y->X)$$
于是二维的Gibbs Sampling采样算法，采样过程如下：

随机初始化$X_0=x_0, Y_0=y_0$
对于t=1,2,…循环采样：
2.1 从条件概率分布$p(y|x_0)$中采样得到$y_1$
2.2 从条件概率分布$p(x|y_1)$中采样得到$x_1$

对于多维的情况，算法也是成立的。例如，一个n维的概率分布$\pi(x_1,x_2,…,x_n)$，可以通过在n个坐标轴上轮换进行采样，得到新的样本集。对于轮换到的任意一个坐标轴$x_i$上的转移，马尔科夫链的状态转移概率为$p(x_i|x_1,x_2,…,x_{i-1},x_{i+1},…,x_n)$。即固定n-1个坐标轴，在某一个坐标轴上移动，同样是满足细致平稳条件。多维的Gibbs Sampling算法采样过程如下：

随机初始化${x_i:i=1,…,n}$
对于t=0,1,2,…循环采样：
2.1 $X_{1}^{t+1} ~p(x_1|x_{2}^{(t)},x_{3}^{(t)},…,x_n^{(t)}$
2.2 $X_{2}^{t+1} ~p(x_2|x_{1}^{(t+1)},x_{3}^{(t)},…,x_n^{(t)}$
2.3 …
2.4 $X_{j}^{t+1} ~p(x_j|x_{1}^{(t+1)},x_{2}^{(t+1)},…,x_{j-1}^{(t+1)},x_{j+1}^{(t)},…,x_n^{(t)}$
2.5 …
2.5 $X_{n}^{t+1} ~p(x_1|x_{1}^{(t+1)},x_{2}^{(t+1)},…,x_{n-1}^{(t+1)}$

同样的，轮换坐标轴不是必须的，可以随机选择某一个坐标轴进行状态转移，只不过常用的Gibbs采样的实现都是基于坐标轴轮换的。