什么是吉布斯采样（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$