1 蒙特卡洛方法
- 1.1 蒙特卡洛方法的作用
- 1.2 非均匀分布采样
- 1.3 分布p(x)不好采样怎么办？
2 什么是吉布斯采样
- 2.1 马尔可夫链
- - 2.1.1 什么是马尔可夫链呢？
  - 2.1.2 为什么我们要引入马尔可夫链？
  - 2.1.3 对给定的分布π\piπ，怎么找到对应的P，使得其为平稳马尔可夫过程
- 2.2 MCMC采样
- 2.3 M-H采样
- 2.4 吉布斯采样（Gibbs）
- - 2.4.1 吉布斯采样原理
  - - 2.4.1.1 二维情况
    - 2.4.1.2 高维情况
  - 2.4.2 吉布斯采样过程
参考资料

1 蒙特卡洛方法

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

1.1 蒙特卡洛方法的作用

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

1.2 非均匀分布采样

上述我们是等间隔采样，其实隐含着假设：x在(a,b)区间是均匀分布的。如果x在(a,b)区间不是均匀分布，而是按照概率p(x)分布的呢？此时利用采样点进行积分估计的结果如下所示：
∫abf(x)dx=∫abf(x)p(x)p(x)dx≈1n∑i=0n−1f(xi)p(xi)\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)}}∫abf(x)dx=∫abp(x)f(x)p(x)dx≈n1i=0∑n−1p(xi)f(xi)
相当于我按照概率p(x)采样n个点，然后根据这n个点对应的f(xi)f(x_i)f(xi)和p(xi)p(x_i)p(xi)值来进行估计。均匀分布相当于一种特例。如下：
p(x)=1b−ap(x)=\frac{1}{b-a}p(x)=b−a1
∫abf(x)dx≈1n∑i=0n−1f(xi)p(xi)=1n∑i=0n−1f(xi)1/(b−a)=b−an∑i=0n−1f(xi)\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)∫abf(x)dx≈n1i=0∑n−1p(xi)f(xi)=n1i=0∑n−11/(b−a)f(xi)=nb−ai=0∑n−1f(xi)

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

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

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

2 什么是吉布斯采样

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

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

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

2.1 马尔可夫链

2.1.1 什么是马尔可夫链呢？

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

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

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

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

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

由细致平稳条件，我们就能得到πP=π\pi P=\piπP=π，推导如下：
∑iπ(i)P(i,j)=∑iπ(j)P(j,i)=π(j)∑iP(j,i)=π(j)\sum_i\pi (i)P(i,j)=\sum_i\pi (j)P(j,i)=\pi (j)\sum_iP(j,i)=\pi (j)i∑π(i)P(i,j)=i∑π(j)P(j,i)=π(j)i∑P(j,i)=π(j)
上式就是πP=π\pi P=\piπP=π 元素展开形式。

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

2.2 MCMC采样

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

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

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

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

MCMC采样过程如下：

已知状态转移矩阵Q，分布π(x)\pi (x)π(x)。设状态转移次数阈值为n1n_1n1，需要样本个数为n2n_2n2。
从任意指定初始状态x0x_0x0,令 t=0t=0t=0.
while t<n1+n2t<n_1+n_2t<n1+n2
a. 从条件分布Q(x∣xt)Q(x|x_t)Q(x∣xt)中采样得到样本x∗x_*x∗
b. 从[0,1]均匀分布采样得到u
c. 若u≤α(xt,x∗)=π(x∗)Q(x∗,xt)u\leq\alpha(x_t,x_*)=\pi (x_*)Q(x_*,x_t)u≤α(xt,x∗)=π(x∗)Q(x∗,xt), 则接受转移，即xt+1=x∗x_{t+1}=x_*xt+1=x∗，t++
d. 否则，不接受转移。t不变。
最后得到的样本集xn1,xn1+1,...,xn1+n2−1x_{n_1},x_{n_1+1},...,x_{n_1+n_2-1}xn1,xn1+1,...,xn1+n2−1就是我们需要的采样结果。

2.3 M-H采样

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

MCMC已经解决了我们的采样问题，但是其样本接受率可能比较低，就是α(xt,x∗)\alpha(x_t,x_*)α(xt,x∗)可能为0.1，0.001，这样采样次数就会比较多，浪费了很多时间。怎样提高接受率呢？
观察细致平稳条件
π(i)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i),其中α(i,j)=π(j)Q(j,i)\pi (i)Q(i,j)\alpha(i,j)=\pi (j)Q(j,i)\alpha(j,i) , 其中\alpha(i,j)=\pi (j)Q(j,i)π(i)Q(i,j)α(i,j)=π(j)Q(j,i)α(j,i),其中α(i,j)=π(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)}
然后两边再同时除以α(j,i)∗α(i,j)\alpha(j,i)*\alpha(i,j)α(j,i)∗α(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}
所以我们得到：

αnew(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是对称的，上式可简化为：
αnew(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采样也有一些问题，比如

当特征比较多时，π(j)Q(j,i)/π(i)Q(i,j)\pi (j)Q(j,i)/\pi (i)Q(i,j)π(j)Q(j,i)/π(i)Q(i,j)的计算比较耗时。

虽然αnew(i,j)\alpha_{new}(i,j)αnew(i,j)的接受率比较大，但还是小于1的，还是有很多辛苦计算的采样点被拒绝掉。

此外，很多时候我们不知道其联合分布，只知道其条件分布，那么这种情况我们就无法对联合分布进行采样。

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

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

2.4.1 吉布斯采样原理

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

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

2.4.1.1 二维情况

假设我们有状态A=(x11,x21)A=(x_1^1, x_2^1)A=(x11,x21), B=(x11,x22)B=(x_1^1, x_2^2)B=(x11,x22), C=(x12,x21)C=(x_1^2,x_2^1)C=(x12,x21), 则有：
π(A)π(x22∣x11)=π(x11,x21)π(x22∣x11)=π(x11)π(x21∣x11)π(x22∣x11)=π(x11)π(x22∣x11)π(x21∣x11)=π(x11,x22)π(x21∣x11)=π(B)π(x21∣x11)\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) π(A)π(x22∣x11)=π(x11,x21)π(x22∣x11)=π(x11)π(x21∣x11)π(x22∣x11)=π(x11)π(x22∣x11)π(x21∣x11)=π(x11,x22)π(x21∣x11)=π(B)π(x21∣x11)
π(A)π(x12∣x21)=π(x11,x21)π(x12∣x21)=π(x21)π(x11∣x21)π(x12∣x21)=π(x21)π(x12∣x21)π(x11∣x21)=π(x12,x21)π(x11∣x21)=π(C)π(x11∣x21)\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) π(A)π(x12∣x21)=π(x11,x21)π(x12∣x21)=π(x21)π(x11∣x21)π(x12∣x21)=π(x21)π(x12∣x21)π(x11∣x21)=π(x12,x21)π(x11∣x21)=π(C)π(x11∣x21)

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

P(A→B)=π(x2B∣x11)，若x1A=x1B=x11P(A \rightarrow B)=\pi(x_2^B|x_1^{1})，若x_1^A= x_1^B= x_1^1P(A→B)=π(x2B∣x11)，若x1A=x1B=x11
P(A→C)=π(x1C∣x21)，若x2A=x2C=x21P(A \rightarrow C)=\pi(x_1^C|x_2^{1})，若x_2^A= x_2^C= x_2^1P(A→C)=π(x1C∣x21)，若x2A=x2C=x21
P(A→D)=0，其他情况P(A \rightarrow D)=0，其他情况P(A→D)=0，其他情况

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

证明如下：
因为满足
x1A=x1B=x11x_1^A= x_1^B= x_1^1x1A=x1B=x11，
所以
P(A→B)=π(x2B∣x11)=π(x22∣x11)P(B→A)=π(x2A∣x11)=π(x21∣x11)P(A \rightarrow B)=\pi(x_2^B|x_1^{1})=\pi(x_2^2|x_1^{1}) \\ P(B \rightarrow A)=\pi(x_2^A|x_1^{1})=\pi(x_2^1|x_1^{1})P(A→B)=π(x2B∣x11)=π(x22∣x11)P(B→A)=π(x2A∣x11)=π(x21∣x11)
所以有：
π(A)P(A→B)=π(A)π(x22∣x11)=π(B)π(x21∣x11)=π(B)P(B→A)\pi(A)P(A\rightarrow B)=\pi (A)\pi(x_2^2|x_1^1) =\pi (B)\pi(x_2^1|x_1^1) =\pi(B)P(B\rightarrow A)π(A)P(A→B)=π(A)π(x22∣x11)=π(B)π(x21∣x11)=π(B)P(B→A)

因为满足
x2A=x2C=x21x_2^A= x_2^C= x_2^1x2A=x2C=x21
所以
P(A→C)=π(x1C∣x21)=π(x12∣x21)P(A \rightarrow C)=\pi(x_1^C|x_2^{1})=\pi(x_1^2|x_2^{1})P(A→C)=π(x1C∣x21)=π(x12∣x21)
P(C→A)=π(x1A∣x21)=π(x11∣x21)P(C \rightarrow A)=\pi(x_1^A|x_2^{1})=\pi(x_1^1|x_2^{1})P(C→A)=π(x1A∣x21)=π(x11∣x21)
所以有：
π(A)P(A→C)=π(A)π(x12∣x21)=π(C)π(x11∣x21)=π(C)π(C→A)\pi(A)P(A\rightarrow C)=\pi (A)\pi(x_1^2|x_2^1)=\pi (C)\pi(x_1^1|x_2^1)=\pi (C)\pi(C \rightarrow A)π(A)P(A→C)=π(A)π(x12∣x21)=π(C)π(x11∣x21)=π(C)π(C→A)

因为既不满足
x1A=x1D=x11x_1^A= x_1^D= x_1^1x1A=x1D=x11
也不满足
x2A=x2D=x21x_2^A= x_2^D= x_2^1x2A=x2D=x21
所以
P(A→D)=0P(A \rightarrow D)=0P(A→D)=0
P(D→A)=0P(D \rightarrow A)=0P(D→A)=0
所以有：
π(A)P(A→D)=0=π(D)π(D→A)\pi(A)P(A\rightarrow D)=0=\pi (D)\pi(D \rightarrow A)π(A)P(A→D)=0=π(D)π(D→A)

综上，得证任意两点E,F, 有
π(E)P(E→F)=π(F)P(F→E)\pi(E)P(E\rightarrow F)=\pi(F)P(F\rightarrow E)π(E)P(E→F)=π(F)P(F→E)

2.4.1.2 高维情况

假设假设A=(x11,x21,...,xn1)A=(x_1^1,x_2^1,...,x_n^1)A=(x11,x21,...,xn1),Bn=(x11,x21,...xn−11,xn2)B_n=(x_1^1,x_2^1,...x_{n-1}^1,x_n^2)Bn=(x11,x21,...xn−11,xn2),我们有：
π(A)π(xn2∣x11,x21,...,xn−11)=π(x11,x21,...,xn1)π(xn2∣x11,x21,...,xn−11)=π(x11,x21,...,xn−11)π(xn1∣x11,x21,...,xn−11)π(xn2∣x11,x21,...,xn−11)=π(x11,x21,...,xn−11)π(xn2∣x11,x21,...,xn−11)π(xn1∣x11,x21,...,xn−11)=π(x11,x21,...xn−11,xn2)π(xn1∣x11,x21,...,xn−11)=π(B)π(xn1∣x11,x21,...,xn−11)\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) π(A)π(xn2∣x11,x21,...,xn−11)=π(x11,x21,...,xn1)π(xn2∣x11,x21,...,xn−11)=π(x11,x21,...,xn−11)π(xn1∣x11,x21,...,xn−11)π(xn2∣x11,x21,...,xn−11)=π(x11,x21,...,xn−11)π(xn2∣x11,x21,...,xn−11)π(xn1∣x11,x21,...,xn−11)=π(x11,x21,...xn−11,xn2)π(xn1∣x11,x21,...,xn−11)=π(B)π(xn1∣x11,x21,...,xn−11)

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

P(A→Bn)=π(xnBn∣x11,x21,...,xn−11)，若xiA=xiBn=xi1,1≤i≤n且i≠n.P(A \rightarrow 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\leq i\leq n且i\neq n.P(A→Bn)=π(xnBn∣x11,x21,...,xn−11)，若xiA=xiBn=xi1,1≤i≤n且i=n.

P(A→Bn−1)=π(xn−1Bn−1∣x11,x21,...,xn−21,xn1)，若xiA=xiBn−1=xi1,1≤i≤n且i≠n−1.P(A \rightarrow 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\leq i\leq n且i\neq n-1.P(A→Bn−1)=π(xn−1Bn−1∣x11,x21,...,xn−21,xn1)，若xiA=xiBn−1=xi1,1≤i≤n且i=n−1.

…

P(A→Bj)=π(xjBj∣x11,x21,...,xj−11,xj+11,...,xn1)，若xiA=xiBj=xi1,1≤i≤n且i≠j.P(A \rightarrow 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\leq i\leq n且i\neq j.P(A→Bj)=π(xjBj∣x11,x21,...,xj−11,xj+11,...,xn1)，若xiA=xiBj=xi1,1≤i≤n且i=j.

…

P(A→B1)=π(x1B1∣x21,x31,...,xn1)，若xiA=xiB1=xi1,1≤i≤n且i≠1.P(A \rightarrow 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\leq i\leq n且i\neq 1.P(A→B1)=π(x1B1∣x21,x31,...,xn1)，若xiA=xiB1=xi1,1≤i≤n且i=1.

P(A→D)=0，其他情况.P(A \rightarrow D)=0，其他情况.P(A→D)=0，其他情况.

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

证明过程同二维情况。

综上，状态转移矩阵P就变成了完全条件概率：
A=(x11,x21,...,xn1)A=(x_1^1,x_2^1,...,x_n^1)A=(x11,x21,...,xn1),
Bj=(x11,x21,...xj−11,xj2,xj+11,...,xn1)B_j=(x_1^1,x_2^1,...x_{j-1}^1,x_{j}^2,x_{j+1}^1,...,x_n^1)Bj=(x11,x21,...xj−11,xj2,xj+11,...,xn1)
P(A→Bj)=π(xjBj∣x11,x21,...,xj−11,xj+11,...,xn1)P(A \rightarrow B_{j})=\pi(x_j^{B_{j}}|x_1^1,x_2^1,...,x_{j-1}^1,x_{j+1}^1,...,x_{n}^1)P(A→Bj)=π(xjBj∣x11,x21,...,xj−11,xj+11,...,xn1)
此时，我们的采样过程就变成了按轴转移的过程了。想更新样本点的第j维，使用第j维的完全条件概率分布π(xjnew∣x1old,x2old,...,xj−1old,xj+1old,...,xnold)\pi(x_j^{new}|x_1^{old},x_2^{old},...,x_{j-1}^{old},x_{j+1}^{old},...,x_{n}^{old})π(xjnew∣x1old,x2old,...,xj−1old,xj+1old,...,xnold),即可采样得到第j维新的值xjnewx_j^{new}xjnew。
所有维度的值都更新一遍，就能得到一个新的样本点(x1new,x2new,...,xnnew)(x_1^{new},x_2^{new},...,x_n^{new})(x1new,x2new,...,xnnew)

2.4.2 吉布斯采样过程

具体采样过程如下：

随机初始化{xi0:i=1,2,...,M}\{x_i^0:i=1,2,...,M\}{xi0:i=1,2,...,M}
for τ=0,1,2,...,n1+n2−1\tau=0,1,2,...,n_1+n_2-1τ=0,1,2,...,n1+n2−1，其中n1n_1n1是采样阈值，n2n_2n2是需要的样本个数。
采样x1τ+1∼p(x1∣x2τ,x3τ,...,xMτ)x_1^{\tau+1} \sim p(x_1|x_2^{\tau},x_3^{\tau},...,x_M^{\tau})x1τ+1∼p(x1∣x2τ,x3τ,...,xMτ)
采样x2τ+1∼p(x2∣x1τ+1,x3τ,...,xMτ)x_2^{\tau+1} \sim p(x_2|x_1^{\tau+1},x_3^{\tau},...,x_M^{\tau})x2τ+1∼p(x2∣x1τ+1,x3τ,...,xMτ)
…
采样xjτ+1∼p(xj∣x1τ+1,x2τ+1,...,xj−1τ+1,xj+1τ...,xMτ)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})xjτ+1∼p(xj∣x1τ+1,x2τ+1,...,xj−1τ+1,xj+1τ...,xMτ)
…
采样xMτ+1∼p(xM∣x1τ+1,x2τ+1,...,xM−1τ+1)x_M^{\tau+1} \sim p(x_M|x_1^{\tau+1},x_2^{\tau+1},...,x_{M-1}^{\tau+1})xMτ+1∼p(xM∣x1τ+1,x2τ+1,...,xM−1τ+1)
以此类推，得到采样点(x1n1,x2n1,...,xMn1),(x1n1+1,x2n1+1,...,xMn1+1),...,(x1n1+n2−1,x2n1+n2−1,...,xMn1+n2−1)(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})(x1n1,x2n1,...,xMn1),(x1n1+1,x2n1+1,...,xMn1+1),...,(x1n1+n2−1,x2n1+n2−1,...,xMn1+n2−1)。

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

参考资料

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

什么是吉布斯采样（Gibbs Sampling）相关推荐

随机采样和随机模拟：吉布斯采样Gibbs Sampling实现高斯分布参数推断
http://blog.csdn.net/pipisorry/article/details/51539739 吉布斯采样的实现问题本文主要说明如何通过吉布斯采样来采样截断多维高斯分布的参数(已知一 ...
随机采样和随机模拟：吉布斯采样Gibbs Sampling
2016年05月12日 00:24:21 阅读数:45658 http://blog.csdn.net/pipisorry/article/details/51373090 吉布斯采样算法详解为什么 ...
三步完成吉布斯采样Gibbs sampling
吉布斯采样的具体执行过程只需要三个步骤,非常非常简单好理解,其它相关的背景知识能帮助加深理解. 一.Preliminaries Monte Carlo methods 它是很宽泛的一类计算方法,依赖重 ...
【ML】线性回归的吉布斯采样(Gibbs Sampling)实现(python)
导航 Bayesian Linear Regression Gibbs Sampling Derving a Gibbs sampler Update for β0\beta_0β0 Update ...
使用MATLAB贝叶斯工具箱(BNT)，进行吉布斯采样(Gibbs Sampling)之前需要做的编译工作...
使用BNT(Bayesian Networks Toolbox)进行推断时,内置了吉布斯采样算法(即gibbs_sampling_inf_engine),但是如果调用这个引擎做推断会报错.报错内容大概 ...
Markov Chain Monte Carlo 和 Gibbs Sampling算法
Welcome To My Blog 一.蒙特卡洛模拟蒙特卡洛模拟(Monte Carlo Simulation)是随机模拟的别名,关于随机模拟的一个重要的问题就是:给定一个概率分布p(x),如何生 ...
吉布斯采样（Gibbs Sampling）
吉布斯采样(Gibbs Sampling) 常用于DBM和DBN,吉布斯采样主要用在像LDA和其它模型参数的推断上. 要完成Gibbs抽样,需要知道条件概率.也就是说,gibbs采样是通过条件 ...
(转载)机器学习知识点(十三)吉布斯采样法(Gibbs Sampling)
3.1 随机模拟随机模拟(或者统计模拟)方法有一个很酷的别名是蒙特卡罗方法(Monte Carlo Simulation).这个方法的发展始于20世纪40年代,和原子弹制造的曼哈顿计划密切相关,当时 ...
Gibbs Sampling（吉布斯采样）
摘要:Gibbs Sampling利用条件概率产生符合分布的样本,用于估计分布的期望,边缘分布:是一种在无法精确计算情况下,用计算机模拟的方法. 什么是Gibbs Sampling Gibbs Sam ...

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

目录