MCMC和Gibbs Sampling

1.随机模拟

       随机模拟又名蒙特卡罗方法，蒙特卡罗方法的源头就是当年用来计算π的著名的的投针实验，由于统计采样的方法成本很高，一直到计算机迅猛发展以后，随机模拟技术才进入实用阶段，对那些确定算法不可行或者不可能解决的问题，蒙特卡罗方法为大家带来希望

       随机模拟技术有一个很重要的问题就是：给定一个概率分布p(x)，如何在计算机中生成它的样本。一般而言，均匀分布U(0,1)的样本是很容易生成的，我们常见的概率分布无论是离散的还是连续的分布都可以基于均匀分布U(0,1)的样本生成，例如正态分布可以基于著名的Box-Muller变换得到，洽谈几个著名的分布，包括指数分布、Gamma分布、t分布、F分布、Beta分布、Dirichlet分布等都可以通过类似的数学变换得到：离散的分布通过均匀分布更加容易生成。

       不过我们并不总是幸运的，当p(x)的分布很复杂，或者p(x)的分布是高维的分布时，样本的生成可能就很困难了，此时就需要一些更加复杂的随机模拟的方法来生成样本，MCMC和Gibbs算法就是最常用的，要了解这两个算法，我们首先对马氏链的平稳分布的性质有基本的认识。

2.马氏链及其平稳分布

      马氏链的数学定义：状态转移的概率只依赖于前一个状态

      在一个例子中：社会学家按照人的经济情况把人分为三类：下层、中层、上层、并且决定一个人收入的重要因素是父母的收入阶层，从父代到子代变化的转移概率是：

转移概率

     使用矩阵方式，转移概率矩阵表示为：

转移概率矩阵

       我们发现当给定不同的初始概率分布时，然后我们计算n代人的概率分布，最后都收敛到相同的概率分布，也就是说，收敛的行为和初始概率分布无关，而是取决于转移概率矩阵P，通过计算，发现P的n次方，发现当n足够大时，P的n次方这个矩阵的每一行都是最后收敛的概率分布，自然的，这个收敛行为并非是这个马氏链所独有的，而是绝大多数马氏链的共同行为，关于马氏链的收敛我们有一个定理：

马氏链收敛定理1

    这个马氏链收敛定理非常重要，所有的MCMC方法都是以这定理为理论基础的，下面给这个定理进行一些说明：

（1）马氏链的状态不要求有限，可以是无穷多个

（2）两个状态i,j连通并非指的是状态i可以一步转移到状态j，而是指的是状态i可以经过有限步转移到状态j，马氏链的任意两个状态连通是指存在一个n，使得矩阵P的n次方中的任何一个元素的数值都大于0

      由马氏链的收敛定理：从初始状态出发，概率分布将收敛到平稳分布，假设经过n步收敛，那么n步以后的转移序列都是平稳分布的样本

3.Markov Chain Monte Carlo

(1)Metropolis算法

      给定概率分布p(x)，我们希望有便捷的方式生成它对应的样本，由于马氏链能收敛到平稳分布，那么：如果我们能构造一个转移矩阵为P的马氏链，使得该马氏链的平稳分布恰好是p(x)，那么我们从任何一个初始状态出发沿着马氏链转移，得到一个转移序列，如果马氏链在第n步收敛了，那么我们取n步以后的转移序列，就得到了平稳分布p(x)的样本

      Metropolis基于这个绝妙的想法，首次提出了基于马氏链的蒙特卡罗方法，即Metropolis算法，并在最早的计算机上编程实现，Metropolis算法是首个普适的采样方法，并启发了一系列的MCMC方法，所以人们视他为随机模拟技术腾飞的起点

(2)MCMC算法

 ①细致平稳条件

     接下来介绍的MCMC算法是Metropolis算法的一个改进变种，即常用的Metropolis-Hastings算法，由前面介绍的可知：马氏链的收敛性质主要由转移矩阵P决定，所以基于马氏链做采样的关键问题是如何构建转移矩阵P，使得平稳分布恰好是我们要的分布p(x)，如何做到这一点呢？我们使用如下的定理：

马氏链收敛定理2

      其实这个定理是显而易见的，因为细致平稳条件的物理意义是：对于任何两个状态i,j，从i转移出去到j而丢失的概率质量，恰好会被j转移回来i丢失的概率质量补充回来，所以状态i上概率质量π(i)是稳定的，从而π(x)是马氏链的平稳分布

②改造马氏链

      假设我们已经有一个转移矩阵Q的马氏链，q(i,j)表示从状态i转移到状态j的概率，显然，通常情况下，p(i)q(i,j)≠p(j)q(j,i)，也就是细致平稳条件不成立，所以p(x)不太可能是这个马氏链的平稳分布。我们可否做一个改造使得细致平稳条件成立呢？

    譬如，我们引入一个α(i,j)，我们希望p(i)q(i,j)α(i,j)=p(j)q(j,i)α(j,i)，那么取什么样的α(i,j)使得这个等式成立呢？最简单的按照对称性，我们取：α(i,j)=p(j)q(j,i)，α(j,i)=p(i)q(i,j)，于是上面这个等式就成立了，所以有：

     于是我们把原来具有转移矩阵Q的一个很普通的马氏链改造为了具有转移矩阵Q'的马氏链，而Q'恰好满足细致平稳条件，由此马氏链Q'的平稳分布就是p(x)!

    在改造Q的过程中引入的α(i,j)称为接受率，物理意义可以理解为在原来的马氏链上，从状态i以q(i,j)的概率跳转到状态j的时候，我们以α(i,j)的概率接受转移，于是得到新的马氏链Q'的转移概率为q(i,j)α(i,j)

马氏链转移和接受概率

③MCMC算法步骤概述

    假设我们已经有一个转移矩阵Q，我们把上面的过程整理一下，就得到了如下采样概率分布p(x)的算法：

MCMC采样算法

1.初始化马氏链初始状态X0=x0,

2.对t=0,1,2...,循环以下过程进行采样

  ●第t个时刻马氏链状态Xt=xt，采样y~q(x|xt)

  ●从均匀分布u~Uniform[0,1]

  ●如果u<α(xt,y)=p(y)q(y,xt)，则接受xt->y,则Xt+1=y

  ●否则不接受转移，即Xt+1=xt

 对这个采样算法的过程说明几点：

      (1)q(x|xt)表示转移概率分布，是在当前状态为xt时，下一个状态的概率分布，取值是可能从xt转移到的状态，而不是概率，对应到转移概率矩阵Q中的元素

     (2)均匀分布随机取样u能落在0到接受率α(xt,y)区间的的概率就是接受率，也就是如果我随机取样取到这个区间，就从当前状态xt转移到下一个状态y，等价于当前状态为xt时按照这个接受率α(xt,y)转移到y，两个说法是一个意思

     上述过程中p(x),q(x|y)说的都是离散的情形，事实上，即使这两个分布是连续的，以上算法仍然是有效的，于是就得到了更一般的连续概率分布p(x)的采样算法，而q(x|y)就是任意一个连续二元概率分布对应的条件分布

      以上的MCMC算法已经能很漂亮的工作了，不过有一个小的问题就是：马氏链Q在转移过程中的接受率α(i,j)可能偏小，这样的采样过程，马氏链容易原地踏步，拒绝大量的跳转，这使得马氏链遍历所有的状态空间可能要花费太长的时间，收敛到平稳分布p(x)的速度太慢，有没有办法提升一些接受率呢？

④改进MCMC算法-Mereopolis-Hastings算法

      假设α(i,j)=0.1，α(j,i)=0.1，此时满足细致平稳条件，于是：p(i)q(i,j)×0.1=p(j)q(j|i)×0.2，将式子 两边扩大五倍，改写为：p(i)q(i,j)×0.5=p(j)q(j|i)×1，看，我们提高了接受率，但是细致平稳条件并没有被打破，这启发我们可以把细致平稳条件p(i)q(i,j)α(i,j)=p(j)q(j,i)α(j,i)两边的α(i,j)，α(j,i)同比例扩大，使得两数中的最大的一个放大到1，这样我们就提高了采样中的跳转接受率，所以我们可以取：

      于是经过对上述MCMC采样算法中接受率的微小改造，我们就得到了最常见的Mereopolis-Hastings算法

      对于分布p(x)，我们构造转移矩阵Q'使其满足细致平稳条件：p(x)Q'(x->y)=p(y)Q'(y->x)，此处x并不要求是一维的，对于高维空间的p(x)，如果满足细致平稳条件：p(X)Q'(X->Y)=p(Y)Q'(Y->X)

Mereopolis-Hastings采样算法

1.初始化马氏链初始状态X0=x0

2.对t=0,1,2...,循环一下过程进行采样

 ●第t个时刻马氏链状态Xt=xt，采样y~q(x|xt)

 ●从均匀分布u~Uniform[0,1]

 ●如果u<α(xt,y)=min{p(y)q(xt|y)/p(xt)q(y|xt),1}，则接受xt->y,则Xt+1=y

 ●否则不接受转移，即Xt+1=xt

那么以上的Mereopolis-Hastings算法一样有效

4.Gibbs Sampling

(1)二维Gibbs Sampling

       对于高维的情形，由于接受率α（通常α<1）的存在，以上Mereopolis-Hastings算法的效率不够高，能否找到一个转移矩阵Q使得接受率α=1呢？

     我们先看看二维的情形，假设有一个概率分布p(x,y)，考察x坐标相同的两个点，A(x1,y1),B(x2,y2)，我们发现，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轴的直线上，如果使用条件分布p(y|x1)作为任何两个点之间的转移概率，那么任何两个点之间的转移满足细致平稳条件。同样的：

     如果我们在y=y1这条直线上任意取两个点A(x1,y1)，C(x2,y1)，也有如下等式：p(A)p(x2|y1)=p(C)p(x1|y1)，也满足细致平稳条件

     那么在沿着坐标轴方向的转移都是满足细致平稳条件的，此时的接受率是1，因为我们不需要两边乘以一个α使得原来的p(i)q(i,j)≠p(j)q(j|i)这个不等式变成等式，因为这个等式在这种情况下是成立的，也可以理解为两边都乘以接受率1，等式成立

平面上马氏链转移矩阵的构造

     于是我们可以构造平面上任意两点之间的转移概率矩阵Q:

     Q(A->B)=p(y2|x1)    如果xA=xB=x1

     Q(A->C)=p(x2|y1)    如果yA=yC=y1

     Q(A->D)=0                其他

     有了如上的转移概率矩阵Q，我们很容易验证对于平面上任意两点X,Y，满足细致平稳条件：p(X)Q(X->Y)=p(Y)Q(Y->X)

    于是这个二维空间上的马氏链将收敛到平稳分布p(x,y)，而这个算法就称为Gibbs Sampling算法

二维Gibbs Sampling算法

1.随机初始化X0=x0,Y0=y0

2.对t=0,1,2...,循环采样

●y(t+1)~p(y|xt)

●x(t+1)~p(x|y(t+1))

       以上采样过程中，马氏链的转移只是轮换沿着坐标轴x轴和y轴做转移，于是得到样本(x0,y0),(x0,y1),(x1,y1),(x1,y2),(x2,y2),....马氏链收敛后，最后得到的样本就是p(x,y)的样本，而收敛之前的阶段称为burn-in period，额外说明一下，教科书上看到的Gibbs Sampling算法大都是沿着坐标轴轮换采样的，但是这是不强求的。

      最一般的情形可以是：在t时刻，可以在x轴和y轴之间随机的选一个坐标轴，然后按照条件概率做转移，马氏链也是一样收敛的，轮换坐标轴只是一种方便的形式

(2)n维Gibbs Sampling

      以上的形式我们很容易推广到高维的情形，如果x1变为,多维情形X1，可以看出推导过程不变，所以细致平稳条件同样是成立的，P(X1,Y1)P(Y2|X1)=P(X1,Y2)P(Y1|X1)，此时转移矩阵Q由条件分布p(Y|X1)定义。上式只是说明一根坐标轴的情形，和二维情形类似，很容易验证对所有坐标轴都有类似的结论，所以n维空间中，对于概率分布p(x1,x2,...,xn)可以定义如下转移矩阵：

     1.如果当前状态是x1,x2,..,xn，马氏链转移的过程中，只能沿着坐标轴做转移。沿着xi这根坐标轴做转移的时候，转移概率由条件概率p(xi|x1,...,xi-1,xi+1,...,xn)定义

     2.其他无法沿着单根坐标轴做的跳转，转移概率都设置为0,

      于是我们把Gibbs Sampling算法从采样二维p(x,y)推广到了采样n维的p(x1,x2,..,xn)

       以上算法收敛后，得到的就是概率分布p(x1,x2,..,xn)的样本，当然这些样本并不独立，但是此处要求的是采样得到的样本符合给定的概率分布，并不要求独立。

      同样的，在以上的算法中，坐标轴轮换采样并不是必须的，可以在坐标轴轮换中引入随机性，这时候转移矩阵Q中任意两个点的转移概率中就会包含坐标轴选择的概率，而在通常的Gibbs Sampling采样算法中，坐标轴轮换是个确定性的过程，也就是在给定时刻t，在一根固定的坐标轴上转移的概率是1

n维Gibbs Sampling算法

1.随机初始化{xi:i=1,2,..,n}

2.对t=0,1,2...,循环采样

5.总结

      本部分的内容主要是介绍几种随机模拟的采样方法，作用是：给定一个概率分布p(x)，然后采样得到这个概率分布的样本，其中发展过程是：

      首先Metropolis基于"如果能构造马氏链的转移概率矩阵P，使得该马氏链的平稳分布恰好是p(x)，那么我们从任何一个初始状态出发沿着马氏链转移，得到一个转移序列，如果马氏链在第n步收敛了，那么我们取第n步以后的即收敛后的序列，就得到了平稳分布p(x)的样本"的思想提出了Metropolis算法，开创了随机模拟技术的起点；

      随后，启发了一系列的MCMC算法，MCMC算法的过程是：马氏链的收敛性质主要由转移矩阵P决定，所以基于马氏链做采样的关键问题是如何构建转移矩阵P，使得平稳分布恰好是我们要的分布p(x)，如何做到这一点呢？首先提出了一个细致平稳条件，根据这个细致平稳条件，我们来构造马氏链的转移矩阵Q，因为不是所有的马氏链都满足这个条件，我们通过引入接受率来保证细致平稳条件成立，即通过重新构造新的转移矩阵Q'使得马氏链满足细致平稳条件，具体实现这一改造马氏链的过程是：首先从任何一个随机初始状态出发，不断的循环采样，采样过程是：每一次状态的跳转取决于接受率：我们从均匀分布U(0,1)取样，如果落在（0，接受率）区间上，那么就按照转移矩阵Q跳转到下一个状态，否则就拒绝跳转，停在当前状态；这样的采样过程存在一个小问题：马氏链Q在转移过程中的接受率α(i,j)可能偏小，这样的采样过程，马氏链容易原地踏步，拒绝大量的跳转，这使得马氏链遍历所有的状态空间可能要花费太长的时间，收敛到平稳分布p(x)的速度太慢，有没有办法提升一些接受率呢？

     此时，我们的新的改进算法Metropolis-Hastings算法出现了，即在MCMC算法的基础上对接受率进行改进，具体是：由细致平稳条件的两边同时扩大同样的倍数，仍然成立，那么我们将等式两边的接受率的最大值增加到1，另外一个接受率增加相同的倍数，在不违背细致平稳条件的情况下，增大了接受率，增加了跳转的几率，即减少了马氏链收敛所需要的时间。问题又来了，由于接受率小于1，Metropolis-Hastings算法在高维情况下的效率不够高，那么能否找到一个接受率为1呢？

      接着又引入一个新的算法是Gibbs Sampling算法，首先是在二维情况下，验证沿着坐标轴方向的转移都是满足细致平稳条件的，此时的接受率是1，在已知概率分布p(x,y)下，基于沿着坐标轴的条件转移概率构造转移矩阵Q，在此转移矩阵下，对于平面上任意两点X,Y，满足细致平稳条件，因此马氏链收敛到平稳分布p(x,y)，采样的过程是从任意一个初始状态出发，沿着状态转移矩阵进行循环采样，每一轮都得到了新的一个状态(x,y)，同样的可延伸到多维情况下。

      到此，按照实际过程介绍了三种随机模拟的采样方法，后面的文章将会介绍文本建模的过程，并将随机模拟采样的方法应用到文本建模中。

参考资料：LDA数学八卦