Gibbs Sampling

algorithm

本文介绍线性同余，平稳马尔可夫，细致平稳过程，Gibbs Sampling

本算法的本质是产生以特定概率产生随机数。

随机模拟与线性同余简介

随机模拟的洋名叫蒙特卡罗方法，目的是通过计算模拟尽可能真实的随机数，对于均匀分布，可以通过线性同余发生器来模拟，其结果与理论值发接近

线性同余简介

同余

$a (mod) m=b(mod)m$，则a与b同余记为$a\equiv b(mod\ m)$,线性同余即为：$ax+c\equiv b\ (mod)\ m$
线性同余序列即为：

$$X_{n+1}=(aX_n+c)\ mod\ m$$
a为乘子，c为增量，m为模

周期

线性同余是有周期的，其周期最大为m，可推理。
定义两个集合S，T，集合S包含从0到m-1所有m个元素，T为一个空集合，那么：
$S=\{0,1,2,\cdots,m-2,m-1\},size(S)=m$
$T=\{\},size(T)=0$
以任意$X_0$为被始值产生$X_1$，它必定在S中，此时将它从S集合中移到T集合中，则：
$S=\{0,1,2,\cdots,m-2,m-1\},size(S)=m-1$
$T=\{X_1\},size(T)=1$
以$X_1$为参数产生第二个随机数时，$X_2$可能属于S，也可能属于T
如果属于S，此时还没产生周期
如果属于T，此时产生周期，周期P=1
假设产生了i-1个都在S中，那么此时S与T的状态为：
$S=\{0,1,2,\cdots,m-2,m-1\},size(S)=m-(i-1)$
$T=\{X_1,X_2,\cdots,X_{i-1}\},size(T)=i-1$
当产生$X_i$时：
如果属于S，此时还没产生周期
如果属于T，此时产生周期，周期P<=i-1
这么推理，一定存在一个周期，且P<=m。
周期与m,a,c有关，在电脑中当然周期是越大越好。更深奥的介绍见：计算机程序设计艺术(The Art of Computer Programming,Volume 2)书太贵！

利用线性同余可以产生常见的分布，如正态分布可以通过Box-Muller变换
设$U_1,U_2\sim Uniform[0,1]$则：
$Z_0=\sqrt{-2\ln{U_1}}\cos(2\pi U_2)$
$Z_1=\sqrt{-2\ln{U_1}}\sin(2\pi U_2)$
则$Z_0,Z_1$独立且服从N（0,1）还有类似指数分布、Gamma 分布、t 分布、F 分布、Beta 分布、Dirichlet 分布等等也都可以用类似的方法生成。书统计模拟中写得具体，但是太贵。
但是当P(x)的形式很复杂时，或者是高维分布时，这个办法就不行了，因此马尔可夫平稳分布就起了作用。

马式链及平其平稳分布

马尔可夫过程也即状态转移，随机过程老师讲得特别详细。举个例子：

这里的1,2,3表示状态下层，中层，上层，转移概率为P。对初始分布$\pi_0=[\pi_0(1),\pi_0(2),\pi_0(3)]$,而下一代分布为：$\pi_1=\pi_0P$,$\pi_n=\pi_{n-1}P=\cdots=\pi_0P^n$。下面列出来一个例子：

再换一个初始概率看看结果

可见，最终$\pi$都不变，而这个性质是由P决定的，称具有这样性质的马氏链为马氏平稳分布。
假设在第n步收敛，那么第n步之后产生的样本就为同分布$\pi(x)$的随机变量。
以上是马尔可夫的基本知识。

Markov Chain Mont Carlo

这里提出的一个想法是：构造出转移矩阵P，使平稳分布恰好是$p(x)$，那么就可以得到一个想要的样本，这个想法得了20世纪很多奖。
细致平稳条件
对分布$\pi(x)和p$满足：

$$\pi(i)P_{ij}=\pi(j)P_{ji}\qquad for\ all\ i,j$$
该分布一定是马氏平稳分布可由下证得：
$$\sum_{i=1}^{\infty}\pi(i)P_{ij}=\sum_{i=1}^{\infty}\pi(j)P_{ji}=\pi(j)\sum_{i=1}^\infty P_{ji}=\pi(j) \Rightarrow \ \pi P=\pi$$
所以$\pi$是平稳分布
假设对一个已有的转移矩阵Q（q(i,j)）表示由i转到j，也可记为q(i|j)，通常
$p(i)q(i,j)\neq p(j)q(j,i)$,也就是细致平稳不成立，这可咋办呢？这里可以对马氏链做一个改造，使细致平稳条件成立，改造如下：
$p(i)q(i,j)\alpha(i,j)=p(j)q(j,i) \alpha(j,i) \tag{*}$
而$\alpha$可以很容易地获得
$\alpha(i,j)=p(j)q(j,i),\qquad \alpha(j,i)=p(i)q(i,j)$
这里，记$Q'(i,j)=q(i,j)\alpha(i,j),Q'(j,i)=q(j,i)\alpha(j,i)$
则*式变为：$p(i)Q'(i,j)=p(j)Q'(j,i)$
这样就人造出来了一个细致平稳分布，那它一定是平稳分布（前面证明过了）。
这里的$\alpha$作用是接收率，当跳转发生时，以$\alpha$的概率接收即可。
由此产生MCMC算法如下图：

该算法有一个小问题，就是如果某些$\alpha$太小，导至接受率太低，那么就可能拒绝大量跳转而原地踏步，解决办法就是两边同时扩大，使大的那一方到1。

如:$\alpha(i,j)=0.1,\alpha(j,i)=0.2$得：

$$p(i)q(i,j)\times0.1=p(j)q(j,i)\times0.2$$
$$同时扩大得：p(i)q(i,j)\times0.5=p(j)q(j,i)\times1$$
$$如此便可得：\alpha(i,j)=\min\{\frac{p(j)q(j,i)}{p(i)q(i,j)},1\}$$
改造后便得Metropolis-Hasting采样算法如下图：

Gibbs Sampling

对于高维情形，由于$\alpha$导致效率不够高，现在是想找到一个Q使$\alpha=1$，先看二维，如图：

对于A，B两点，可见：$p(x1,y1)p(y2|x1)=p(x1)p(y1|x1)p(y2|x1)$，$p(x1,y2)p(y1|x1)=p(x1)p(y2|x1)p(y1|x1)$，所以可得：

$$p(x1,y1)p(y2|x1)=p(x1,y2)p(y1|x1)$$
$$即：P(A)p(y2|x1)=p(B)p(y1|x1)$$
可见，在x=x1这条线上任意两个点之间满足细至平稳，同样，y也类似。如此便可构造平面上的转移矩阵Q：
$$\begin{aligned} &Q(A\rightarrow B)=p(y_B|x_1)&&如果&x_A=x_B=x_1\\ &Q(A\rightarrow C)=p(x_c|y_1)&&如果&y_A=y_C=y_1\\ &Q(A\rightarrow D)=0 &&&其他 \end{aligned}$$

由上可得平面两点满足细致平稳条件：

$$P(X)Q(X\rightarrow Y)=p(Y)Q(Y\rightarrow X)$$
对于平面上这个二维点进行平稳分布采样，称为Gibbs Sampling,如下图

如图所示，马氏链转移只沿坐标轴转移，收敛后得到的样本就是p(x,y)样本。其实采样并不需要按坐标轴轮换采样。

还可推广到多维情形，如对n维 概率分布$p(x_1,x_2,\cdots,x_n)$,1.如果当前状态为$(x_1,x_2,\cdots,x_n)$，转移时由条件概率$p(x_i|x_1,\cdots,x_{i-1},X_{i+1},\cdots,x_n)$。2.其他无法沿着单根坐标轴进行的跳转，转移概率都设为0。
算法如下图：