随机模拟（MCMC）

http://cos.name/2013/01/lda-math-mcmc-and-gibbs-sampling/

http://blog.csdn.net/lin360580306/article/details/51240398

http://blog.csdn.net/pipisorry/article/details/51373090

随机模拟(或者统计模拟)方法有一个很酷的别名是蒙特卡罗方法(Monte Carlo Simulation)。这个方法的发展始于20世纪40年代，和原子弹制造的曼哈顿计划密切相关，当时的几个大牛，包括乌拉姆、冯.诺依曼、费米、费曼、Nicholas Metropolis， 在美国洛斯阿拉莫斯国家实验室研究裂变物质的中子连锁反应的时候，开始使用统计模拟的方法,并在最早的计算机上进行编程实现。

随机模拟与计算机

现代的统计模拟方法最早由数学家乌拉姆提出，被Metropolis命名为蒙特卡罗方法，蒙特卡罗是著名的赌场，赌博总是和统计密切关联的，所以这个命名风趣而贴切，很快被大家广泛接受。被不过据说费米之前就已经在实验中使用了，但是没有发表。说起蒙特卡罗方法的源头，可以追溯到18世纪，布丰当年用于计算

蒙特卡罗方法

统计模拟中有一个重要的问题就是给定一个概率分布Uniform(0,1) 的理论计算结果非常接近。这样的伪随机序列就有比较好的统计性质，可以被当成真实的随机数使用。

生成一个概率分布的样本

而我们常见的概率分布，无论是连续的还是离散的分布，都可以基于Uniform(0,1) 的样本生成。例如正态分布可以通过著名的 Box-Muller 变换得到

[Box-Muller 变换]  如果随机变量 U1,U2∼Uniform[0,1]，

Z0=−2ln⁡U1cos(2πU2)Z1=−2ln⁡U1sin(2πU2)

则 Z0,Z1 独立且服从标准正态分布。

 

其它几个著名的连续分布，包括指数分布、Gamma 分布、t 分布、F 分布、Beta 分布、Dirichlet 分布等等,也都可以通过类似的数学变换得到；离散的分布通过均匀分布更加容易生成。更多的统计分布如何通过均匀分布的变换生成出来，大家可以参考统计计算的书，其中 Sheldon M. Ross 的《统计模拟》是写得非常通俗易懂的一本。

不过我们并不是总是这么幸运的，当p(x) 是个高维的分布的时候，样本的生成就可能很困难了。 譬如有如下的情况

• p~(x) 我们是可以计算的，但是底下的积分式无法显式计算。
• p(x) 是高维的，这种情形就更加明显。

此时就需要使用一些更加复杂的随机模拟的方法来生成样本。而本节中将要重点介绍的 MCMC(Markov Chain Monte Carlo) 和 Gibbs Sampling算法就是最常用的一种，这两个方法在现代贝叶斯分析中被广泛使用。要了解这两个算法，我们首先要对马氏链的平稳分布的性质有基本的认识。

3.2 马氏链及其平稳分布

马氏链的数学定义很简单

P(Xt+1=x|Xt,Xt−1,⋯)=P(Xt+1=x|Xt)

也就是状态转移的概率只依赖于前一个状态。

 

我们先来看马氏链的一个具体的例子。社会学家经常把人按其经济状况分成3类：下层(lower-class)、中层(middle-class)、上层(upper-class)，我们用1,2,3 分别代表这三个阶层。社会学家们发现决定一个人的收入阶层的最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别，那么他的孩子属于下层收入的概率是 0.65, 属于中层收入的概率是 0.28, 属于上层收入的概率是 0.07。事实上，从父代到子代，收入阶层的变化的转移概率如下

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

P=[0.650.280.070.150.670.180.120.360.52]

 

假设当前这一代人处在下层、中层、上层的人的比例是概率分布向量 πn=πn−1P=π0Pn。

假设初始概率分布为n代人的分布状况如下

我们发现从第7代人开始，这个分布就稳定不变了，这个是偶然的吗？我们换一个初始概率分布n代人的分布状况如下

我们发现，到第9代人的时候, 分布又收敛了。最为奇特的是，两次给定不同的初始概率分布，最终都收敛到概率分布 Pn

 

P20=P21=⋯=P100=⋯=[0.2860.4890.2250.2860.4890.2250.2860.4890.225]

 

我们发现，当 π=[0.286,0.489,0.225] 这个概率分布。自然的，这个收敛现象并非是我们这个马氏链独有的，而是绝大多数马氏链的共同行为，关于马氏链的收敛我们有如下漂亮的定理：

马氏链定理： 如果一个非周期马氏链具有转移概率矩阵limn→∞Pijn=π(j), 我们有

1. limn→∞Pn=[π(1)π(2)⋯π(j)⋯π(1)π(2)⋯π(j)⋯⋯⋯⋯⋯⋯π(1)π(2)⋯π(j)⋯⋯⋯⋯⋯⋯]

1.  π(j)=∑i=0∞π(i)Pij
2.  πP=π 的唯一非负解

其中,

π=[π(1),π(2),⋯,π(j),⋯],∑i=0∞πi=1

π称为马氏链的平稳分布。

 

这个马氏链的收敛定理非常重要，所有的 MCMC(Markov Chain Monte Carlo) 方法都是以这个定理作为理论基础的。 定理的证明相对复杂，一般的随机过程课本中也不给证明，所以我们就不用纠结它的证明了，直接用这个定理的结论就好了。我们对这个定理的内容做一些解释说明：

1. 该定理中马氏链的状态不要求有限，可以是有无穷多个的；
2. 定理中的“非周期“这个概念我们不打算解释了，因为我们遇到的绝大多数马氏链都是非周期的；
3. 两个状态Pn 中的任何一个元素的数值都大于零。
4. 我们用 limn→∞Pijn=π(j)存在，很容易证明以上定理的第二个结论。由于
   
   P(Xn+1=j)=∑i=0∞P(Xn=i)P(Xn+1=j|Xn=i)=∑i=0∞P(Xn=i)Pij
   π(j)=∑i=0∞π(i)Pij

从初始概率分布 πi, 则有

X0∼π0(x)Xi∼πi(x),πi(x)=πi−1(x)P=π0(x)Pn

由马氏链收敛的定理, 概率分布n步的时候马氏链收敛，则有

X0∼π0(x)X1∼π1(x)⋯Xn∼πn(x)=π(x)Xn+1∼π(x)Xn+2∼π(x)⋯

所以 π(x) 的样本。

 

3.3 Markov Chain Monte Carlo

对于给定的概率分布xn,xn+1⋯。

这个绝妙的想法在1953年被 Metropolis想到了，为了研究粒子系统的平稳性质， Metropolis 考虑了物理学中常见的波尔兹曼分布的采样问题，首次提出了基于马氏链的蒙特卡罗方法，即Metropolis算法，并在最早的计算机上编程实现。Metropolis 算法是首个普适的采样方法，并启发了一系列 MCMC方法，所以人们把它视为随机模拟技术腾飞的起点。 Metropolis的这篇论文被收录在《统计学中的重大突破》中， Metropolis算法也被遴选为二十世纪的十个最重要的算法之一。

我们接下来介绍的MCMC 算法是 Metropolis 算法的一个改进变种，即常用的 Metropolis-Hastings 算法。由上一节的例子和定理我们看到了，马氏链的收敛性质主要由转移矩阵p(x)。如何能做到这一点呢？我们主要使用如下的定理。

定理：[细致平稳条件] 如果非周期马氏链的转移矩阵π(x) 满足

(1)π(i)Pij=π(j)Pjifor alli,j

则 π(x) 是马氏链的平稳分布，上式被称为细致平稳条件(detailed balance condition)。

 

其实这个定理是显而易见的，因为细致平稳条件的物理含义就是对于任何两个状态π(x)是马氏链的平稳分布。数学上的证明也很简单，由细致平稳条件可得

∑i=1∞π(i)Pij=∑i=1∞π(j)Pji=π(j)∑i=1∞Pji=π(j)⇒πP=π

由于π是平稳分布。

 

假设我们已经有一个转移矩阵为q(i→j)), 显然，通常情况下

p(i)q(i,j)≠p(j)q(j,i)

也就是细致平稳条件不成立，所以 α(i,j), 我们希望

(2)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)

于是(*)式就成立了。所以有

(3)p(i)q(i,j)α(i,j)⏟Q′(i,j)=p(j)q(j,i)α(j,i)⏟Q′(j,i)(∗∗)

于是我们把原来具有转移矩阵p(x)！

 

在改造 q(i,j)α(i,j)。

马氏链转移和接受概率

假设我们已经有一个转移矩阵Q(对应元素为p(x)的算法。

上述过程中 q(x|y) 就是任意一个连续二元概率分布对应的条件分布。

以上的 MCMC 采样算法已经能很漂亮的工作了，不过它有一个小的问题：马氏链p(x)的速度太慢。有没有办法提升一些接受率呢?

假设 α(i,j)=0.1,α(j,i)=0.2, 此时满足细致平稳条件，于是

p(i)q(i,j)×0.1=p(j)q(j,i)×0.2

上式两边扩大5倍，我们改写为

p(i)q(i,j)×0.5=p(j)q(j,i)×1

看，我们提高了接受率，而细致平稳条件并没有打破！这启发我们可以把细致平稳条件(**) 式中的α(i,j),α(j,i) 同比例放大，使得两数中最大的一个放大到1，这样我们就提高了采样中的跳转接受率。所以我们可以取

α(i,j)=min{p(j)q(j,i)p(i)q(i,j),1}

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

 

对于分布 Q′ 使其满足细致平稳条件

p(x)Q′(x→y)=p(y)Q′(y→x)

此处 p(x)，如果满足细致平稳条件

p(x)Q′(x→y)=p(y)Q′(y→x)

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

 

3.2 Gibbs Sampling

对于高维的情形，由于接受率 A(x1,y1),B(x1,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)

所以得到

(4)p(x1,y1)p(y2|x1)=p(x1,y2)p(y1|x1)(∗∗∗)

即

p(A)p(y2|x1)=p(B)p(y1|x1)

基于以上等式，我们发现，在 A(x1,y1),C(x2,y1),也有如下等式

p(A)p(x2|y1)=p(C)p(x1|y1).

 

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

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

Q(A→B)=p(yB|x1)如果xA=xB=x1Q(A→C)=p(xC|y1)如果yA=yC=y1Q(A→D)=0其它

 

有了如上的转移矩阵 Q, 我们很容易验证对平面上任意两点 X,Y, 满足细致平稳条件

p(X)Q(X→Y)=p(Y)Q(Y→X)

于是这个二维空间上的马氏链将收敛到平稳分布 p(x,y)。而这个算法就称为 Gibbs Sampling 算法,是 Stuart Geman 和Donald Geman 这两兄弟于1984年提出来的，之所以叫做Gibbs Sampling 是因为他们研究了Gibbs random field, 这个算法在现代贝叶斯分析中占据重要位置。

 

Gibbs Sampling 算法中的马氏链转移

以上采样过程中，如图所示，马氏链的转移只是轮换的沿着坐标轴 y轴之间随机的选一个坐标轴，然后按条件概率做转移，马氏链也是一样收敛的。轮换两个坐标轴只是一种方便的形式。

以上的过程我们很容易推广到高维的情形，对于(***) 式，如果x1，可以看出推导过程不变，所以细致平稳条件同样是成立的

(5)p(x1,y1)p(y2|x1)=p(x1,y2)p(y1|x1)

此时转移矩阵 Q 由条件分布 p(x1,x2,⋯,xn) 可以如下定义转移矩阵

 

1. 如果当前状态为p(xi|x1,⋯,xi−1,xi+1,⋯,xn) 定义；
2.  其它无法沿着单根坐标轴进行的跳转，转移概率都设置为 0。

于是我们可以把Gibbs Smapling 算法从采样二维的 p(x1,x2,⋯,xn)

以上算法收敛后，得到的就是概率分布t，在一根固定的坐标轴上转移的概率是1。