LDA-math-MCMC 和 Gibbs Sampling

关键词：Gibbs 抽样; MCMC; Monte-Carlo; 统计模拟

https://site.douban.com/182577/widget/notes/10567181/note/292072927/
https://cosx.org/2013/01/lda-math-mcmc-and-gibbs-sampling

3.1 随机模拟

蒙特卡罗方法
生成一个概率分布的样本
[Box-Muller 变换]
如果随机变量 U1,U2U1,U2 独立且U1,U2∼Uniform[0,1]U1,U2∼Uniform[0,1]，
Z0=√−2lnU1cos(2πU2)Z1=√−2lnU1sin(2πU2)Z0=−2ln⁡U1cos(2πU2)Z1=−2ln⁡U1sin(2πU2)
则 Z0,Z1Z0,Z1 独立且服从标准正态分布。

其它几个著名的连续分布，包括指数分布、Gamma 分布、t 分布、F 分布、Beta 分布、Dirichlet 分布等等, 也都可以通过类似的数学变换得到；离散的分布通过均匀分布更加容易生成。更多的统计分布如何通过均匀分布的变换生成出来，大家可以参考统计计算的书，其中 Sheldon M. Ross 的《统计模拟》是写得非常通俗易懂的一本。
不过我们并不是总是这么幸运的，当p(x)p(x)的形式很复杂，或者 p(x)p(x) 是个高维的分布的时候，样本的生成就可能很困难了。 譬如有如下的情况

• p(x)=~p(x)∫~p(x)dxp(x)=p~(x)∫p~(x)dx, 而 ~p(x)p~(x) 我们是可以计算的，但是底下的积分式无法显式计算。
• p(x,y)p(x,y) 是一个二维的分布函数，这个函数本身计算很困难，但是条件分布 p(x|y),p(y|x)p(x|y),p(y|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)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⎤⎥⎦P=[0.650.280.070.150.670.180.120.360.52]
假设当前这一代人处在下层、中层、上层的人的比例是概率分布向量 π0=[π0(1),π0(2),π0(3)]π0=[π0(1),π0(2),π0(3)]，那么他们的子女的分布比例将是 π1=π0Pπ1=π0P, 他们的孙子代的分布比例将是 π2=π1P=π0P2π2=π1P=π0P2, ……, 第nn代子孙的收入分布比例将是 πn=πn−1P=π0Pnπn=πn−1P=π0Pn。
假设初始概率分布为π0=[0.21,0.68,0.11]π0=[0.21,0.68,0.11]，则我们可以计算前nn代人的分布状况如下

我们发现从第 7 代人开始，这个分布就稳定不变了，这个是偶然的吗？我们换一个初始概率分布π0=[0.75,0.15,0.1]π0=[0.75,0.15,0.1] 试试看，继续计算前nn代人的分布状况如下

我们发现，到第 9 代人的时候, 分布又收敛了。最为奇特的是，两次给定不同的初始概率分布，最终都收敛到概率分布 π=[0.286,0.489,0.225]π=[0.286,0.489,0.225]，也就是说收敛的行为和初始概率分布 π0π0 无关。这说明这个收敛行为主要是由概率转移矩阵 $P$ 决定的。
我们发现，当 nn 足够大的时候，这个PnPn矩阵的每一行都是稳定地收敛到π=[0.286,0.489,0.225]π=[0.286,0.489,0.225] 这个概率分布。自然的，这个收敛现象并非是我们这个马氏链独有的，而是绝大多数马氏链的共同行为，关于马氏链的收敛我们有如下漂亮的定理：
马氏链定理：

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

1. 该定理中马氏链的状态不要求有限，可以是有无穷多个的；
2. 定理中的 “非周期” 这个概念我们不打算解释了，因为我们遇到的绝大多数马氏链都是非周期的；
3. 两个状态i,ji,j是连通并非指ii 可以直接一步转移到jj(Pij>0Pij>0)，而是指 ii 可以通过有限的nn步转移到达jj(Pnij>0Pijn>0)。马氏链的任何两个状态是连通的含义是指存在一个nn，使得矩阵PnPn 中的任何一个元素的数值都大于零。
4. 我们用 XiXi 表示在马氏链上跳转第ii步后所处的状态，如果limn→∞Pnij=π(j)limn→∞Pijn=π(j)存在，很容易证明以上定理的第二个结论。由于

P(Xn+1=j)=∞∑i=0P(Xn=i)P(Xn+1=j|Xn=i)=∞∑i=0P(Xn=i)PijP(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π(j)=∑i=0∞π(i)Pij
从初始概率分布π0π0出发，我们在马氏链上做状态转移，记XiXi的概率分布为πiπi，则有
X0∼π0(x)Xi∼πi(x),πi(x)=πi−1(x)P=π0(x)PnX0∼π0(x)Xi∼πi(x),πi(x)=πi−1(x)P=π0(x)Pn
由马氏链收敛的定理, 概率分布πi(x)πi(x)将收敛到平稳分布 π(x)π(x)。假设到第nn步的时候马氏链收敛，则有
X0∼π0(x)X1∼π1(x)⋯Xn∼πn(x)=π(x)Xn+1∼π(x)Xn+2∼π(x)⋯X0∼π0(x)X1∼π1(x)⋯Xn∼πn(x)=π(x)Xn+1∼π(x)Xn+2∼π(x)⋯
所以 Xn,Xn+1,Xn+2,⋯∼π(x)Xn,Xn+1,Xn+2,⋯∼π(x) 都是同分布的随机变量，当然他们并不独立。如果我们从一个具体的初始状态 x0x0 开始, 沿着马氏链按照概率转移矩阵做跳转，那么我们得到一个转移序列 x0,x1,x2,⋯xn,xn+1⋯,x0,x1,x2,⋯xn,xn+1⋯, 由于马氏链的收敛行为， xn,xn+1,⋯xn,xn+1,⋯ 都将是平稳分布 π(x)π(x) 的样本。

3.3 Markov Chain Monte Carlo

对于给定的概率分布p(x)p(x), 我们希望能有便捷的方式生成它对应的样本。由于马氏链能收敛到平稳分布， 于是一个很的漂亮想法是：如果我们能构造一个转移矩阵为 $P$ 的马氏链，使得该马氏链的平稳分布恰好是p(x)p(x)， 那么我们从任何一个初始状态x0x0出发沿着马氏链转移, 得到一个转移序列x0,x1,x2,⋯xn,xn+1⋯,x0,x1,x2,⋯xn,xn+1⋯,， 如果马氏链在第nn步已经收敛了，于是我们就得到了π(x)π(x)的样本xn,xn+1⋯xn,xn+1⋯。
这个绝妙的想法在 1953 年被 Metropolis 想到了，为了研究粒子系统的平稳性质， Metropolis 考虑了物理学中常见的波尔兹曼分布的采样问题，首次提出了基于马氏链的蒙特卡罗方法，即 Metropolis 算法，并在最早的计算机上编程实现。Metropolis 算法是首个普适的采样方法，并启发了一系列 MCMC 方法，所以人们把它视为随机模拟技术腾飞的起点。 Metropolis 的这篇论文被收录在《统计学中的重大突破》中， Metropolis 算法也被遴选为二十世纪的十个最重要的算法之一。
我们接下来介绍的 MCMC 算法是 Metropolis 算法的一个改进变种，即常用的 Metropolis-Hastings 算法。由上一节的例子和定理我们看到了，马氏链的收敛性质主要由转移矩阵PP决定, 所以基于马氏链做采样的关键问题是如何构造转移矩阵PP, 使得平稳分布恰好是我们要的分布p(x)p(x)。如何能做到这一点呢？我们主要使用如下的定理。
定理：[细致平稳条件] 如果非周期马氏链的转移矩阵PP和分布π(x)π(x)满足
π(i)Pij=π(j)Pjifor alli,jπ(i)Pij=π(j)Pjifor alli,j
则π(x)π(x)是马氏链的平稳分布，上式被称为细致平稳条件 (detailed balance condition)。
其实这个定理是显而易见的，因为细致平稳条件的物理含义就是对于任何两个状态i,ji,j，从 ii 转移出去到jj而丢失的概率质量，恰好会被从jj转移回ii的概率质量补充回来，所以状态 $i$ 上的概率质量π(i)π(i)是稳定的，从而π(x)π(x)是马氏链的平稳分布。
马氏链转移和接受概率
假设我们已经有一个转移矩阵 Q(对应元素为q(i,j)q(i,j)), 把以上的过程整理一下，我们就得到了如下的用于采样概率分布p(x)p(x)的算法。


对于分布 p(x)p(x), 我们构造转移矩阵 Q′Q′ 使其满足细致平稳条件
p(x)Q′(x→y)=p(y)Q′(y→x)p(x)Q′(x→y)=p(y)Q′(y→x)
此处xx并不要求是一维的，对于高维空间的 p(x)p(x)，如果满足细致平稳条件
p(x)Q′(x→y)=p(y)Q′(y→x)p(x)Q′(x→y)=p(y)Q′(y→x)
那么以上的 Metropolis-Hastings 算法一样有效。

3.4 Gibbs Sampling

对于高维的情形，由于接受率αα的存在 (通常 α<1α<1)，以上 Metropolis-Hastings 算法的效率不够高。能否找到一个转移矩阵 Q 使得接受率 α=1α=1 呢？我们先看看二维的情形，假设有一个概率分布 p(x,y)p(x,y)， 考察 $x$ 坐标相同的两个点A(x1,y1),B(x1,y2)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)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(x1,y1)p(y2|x1)=p(x1,y2)p(y1|x1)(∗∗∗)
即
p(A)p(y2|x1)=p(B)p(y1|x1)p(A)p(y2|x1)=p(B)p(y1|x1)
基于以上等式，我们发现，在 x=x1x=x1 这条平行于 yy轴的直线上，如果使用条件分布 p(y|x1)p(y|x1)做为任何两个点之间的转移概率，那么任何两个点之间的转移满足细致平稳条件。同样的，如果我们在y=y1y=y1 这条直线上任意取两个点 A(x1,y1),C(x2,y1)A(x1,y1),C(x2,y1)，也有如下等式
p(A)p(x2|y1)=p(C)p(x1|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(A→B)=p(yB|x1)如果xA=xB=x1Q(A→C)=p(xC|y1)如果yA=yC=y1Q(A→D)=0其它
有了如上的转移矩阵 Q，我们很容易验证对平面上任意两点 X,YX,Y， 满足细致平稳条件
p(X)Q(X→Y)=p(Y)Q(Y→X)p(X)Q(X→Y)=p(Y)Q(Y→X)
于是这个二维空间上的马氏链将收敛到平稳分布 p(x,y)p(x,y)。而这个算法就称为 Gibbs Sampling 算法，是 Stuart Geman 和 Donald Geman 这两兄弟于 1984 年提出来的，之所以叫做 Gibbs Sampling 是因为他们研究了 Gibbs random field，这个算法在现代贝叶斯分析中占据重要位置。


Gibbs Sampling 算法中的马氏链转移
以上采样过程中，如图所示，马氏链的转移只是轮换的沿着坐标轴xx轴和yy轴做转移，于是得到样本(x0,y0),(x0,y1),(x1,y1),(x1,y2),(x2,y2),⋯(x0,y0),(x0,y1),(x1,y1),(x1,y2),(x2,y2),⋯马氏链收敛后，最终得到的样本就是 p(x,y)p(x,y)的样本，而收敛之前的阶段称为 burn-in period。额外说明一下，我们看到教科书上的 Gibbs Sampling 算法大都是坐标轴轮换采样的，但是这其实是不强制要求的。最一般的情形可以是，在tt时刻，可以在xx轴和yy轴之间随机的选一个坐标轴，然后按条件概率做转移，马氏链也是一样收敛的。轮换两个坐标轴只是一种方便的形式。
以上的过程我们很容易推广到高维的情形，对于 (*) 式，如果x1x1变为多维情形x1x1，可以看出推导过程不变，所以细致平稳条件同样是成立的
p(x1,y1)p(y2|x1)=p(x1,y2)p(y1|x1)p(x1,y1)p(y2|x1)=p(x1,y2)p(y1|x1)
此时转移矩阵 Q 由条件分布p(y|x1)p(y|x1)定义。上式只是说明了一根坐标轴的情形，和二维情形类似，很容易验证对所有坐标轴都有类似的结论。所以 $n$ 维空间中对于概率分布 p(x1,x2,⋯,xn)p(x1,x2,⋯,xn)可以如下定义转移矩阵
1. 如果当前状态为(x1,x2,⋯,xn)(x1,x2,⋯,xn)，马氏链转移的过程中，只能沿着坐标轴做转移。沿着xixi这根坐标轴做转移的时候，转移概率由条件概率p(xi|x1,⋯,xi−1,xi+1,⋯,xn)p(xi|x1,⋯,xi−1,xi+1,⋯,xn)定义； 2. 其它无法沿着单根坐标轴进行的跳转，转移概率都设置为 0。
于是我们可以把 Gibbs Smapling 算法从采样二维的p(x,y)p(x,y)推广到采样nn维的p(x1,x2,⋯,xn)p(x1,x2,⋯,xn)

以上算法收敛后，得到的就是概率分布p(x1,x2,⋯,xn)p(x1,x2,⋯,xn)的样本，当然这些样本并不独立，但是我们此处要求的是采样得到的样本符合给定的概率分布，并不要求独立。同样的，在以上算法中，坐标轴轮换采样不是必须的，可以在坐标轴轮换中引入随机性，这时候转移矩阵QQ中任何两个点的转移概率中就会包含坐标轴选择的概率，而在通常的 Gibbs Sampling 算法中，坐标轴轮换是一个确定性的过程，也就是在给定时刻tt，在一根固定的坐标轴上转移的概率是 1。

MCMC（马尔科夫链蒙特卡洛）

线性同余发生器可以生成伪随机数
线性同余方法（LCG）
它是根据递归公式：

其中是产生器设定的常数。
LCG的周期最大为，但大部分情况都会少于M。要令LCG达到最大周期，应符合以下条件：

1. 互质；
2. 的所有质因子的积能整除；
3. 若是4的倍数，也是；
4. 都比小；

5. 是正整数。

   线性同余算法有m 、a 、c 和X0 4个参数,通过置Xn + 1 ≡aXn + c (mod m) ,求得随机数序列< Xn > , 这个序列称作线性同余序列。m、a 、c 和X0 分别称做模数、乘数、增量和初始值。线性同余方法速度快,如果对乘数和模数进行适当的选择,可以满足用于评价一个随机数产生器的3 种准则:
   1.这个函数应该是一个完整周期的产生函数。也就是说,这个函数应该在重复之前产生出0 到m之间的所有数;
   2.产生的序列应该看起来是随机的;
   3.这个函数应该用32bit 算术高效实现。
   产生随机种子的方法很多，目前用得比较多的是使用系统时间为种子。我觉得这种方法也不妥当。假如我批量执行程序，程序执行的时间是几个ms，那么几个相邻程序的种子就是一样的，产生的结果因此也是一样的。（因为系统时间是按照秒来计算的，一秒内执行多少次，产生的随机种子就有多少相同的。）
   (随机数是否是真随机数？？？)
   生成一个概率分布的样本
   而我们常见的概率分布，无论是连续的还是离散的分布，都可以基于Uniform(0,1)的样本生成。例如正态分布可以通过著名的 Box-Muller 变换得到
   [Box-Muller 变换]
   马氏链及其平稳分布
   马氏链的数学定义很简单
   P(Xt+1=x|Xt,Xt−1,⋯)=P(Xt+1=x|Xt)
   也就是状态转移的概率只依赖于前一个状态。
   马氏链定理： 如果一个非周期马氏链具有转移概率矩阵P,且它的任何两个状态是连通的，那么 limn→∞Pnij 存在且与i无关，记 limn→∞Pnij=π(j), 我们有

6. limn→∞Pn=⎡⎣⎢⎢⎢⎢⎢π(1)π(1)⋯π(1)⋯π(2)π(2)⋯π(2)⋯⋯⋯⋯⋯⋯π(j)π(j)⋯π(j)⋯⋯⋯⋯⋯⋯⎤⎦⎥⎥⎥⎥⎥

7. π(j)=∑i=0∞π(i)Pij
8. π 是方程 πP=π 的唯一非负解

其中,
π=[π(1),π(2),⋯,π(j),⋯],∑i=0∞πi=1
π称为马氏链的平稳分布。

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

1. 该定理中马氏链的状态不要求有限，可以是有无穷多个的；
2. 定理中的“非周期“这个概念我们不打算解释了，因为我们遇到的绝大多数马氏链都是非周期的；
3. 两个状态i,j是连通并非指i 可以直接一步转移到j(Pij>0),而是指 i 可以通过有限的n步转移到达j(Pnij>0)。马氏链的任何两个状态是连通的含义是指存在一个n, 使得矩阵Pn 中的任何一个元素的数值都大于零。
4. 我们用 Xi 表示在马氏链上跳转第i步后所处的状态，如果 limn→∞Pnij=π(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
从初始概率分布 π0 出发，我们在马氏链上做状态转移，记Xi的概率分布为πi, 则有
X0Xi∼π0(x)∼πi(x),πi(x)=πi−1(x)P=π0(x)Pn
由马氏链收敛的定理, 概率分布πi(x)将收敛到平稳分布 π(x)。假设到第n步的时候马氏链收敛，则有
X0X1XnXn+1Xn+2∼π0(x)∼π1(x)⋯∼πn(x)=π(x)∼π(x)∼π(x)⋯
所以 Xn,Xn+1,Xn+2,⋯∼π(x) 都是同分布的随机变量，当然他们并不独立。如果我们从一个具体的初始状态 x0 开始,沿着马氏链按照概率转移矩阵做跳转，那么我们得到一个转移序列 x0,x1,x2,⋯xn,xn+1⋯, 由于马氏链的收敛行为， xn,xn+1,⋯都将是平稳分布 π(x) 的样本。
(1)贝叶斯推断的计算问题

在上节我们看到，贝叶斯统计学是利用后验分布对θ进行推断。这种推断的计算很多情况下要用积分计算来完成。比如，我们要计算θ的函数g(θ)的期望：

E(g(θ∣x))=∫g(θ)fθ∣x(θ∣x)dθ

其中函数f表示后验分布。当g(θ)=θ时，得到的就是关于θ的点估计。

但是对很多贝叶斯推断问题来说，有时候后验分布过于复杂，使得积分没有显示结果，数值方法也很难应用；有时候需要计算多重积分（比如后验分布是多元分布时）。这些都会带来计算上的很大困难。这也是在很长的时期内，贝叶斯统计得不到快速发展的一个原因。1990年代MCMC（Markov Chain Monte Carlo ，马尔科夫链蒙特卡洛）计算方法引入到贝叶斯统计学之后，一举解决了这个计算的难题。可以说，近年来贝叶斯统计的蓬勃发展，特别是在各个学科的广泛应用和MCMC方法的使用有着极其密切的关系。

（2）蒙特卡洛方法（Monte Carlo）

蒙特卡洛方法是一种随机模拟方法，随机模拟的思想由来已久（参见下面的蒲丰投针的例子），但是由于难于取得随机数，随机模拟的方法一直发展缓慢。而蒙特卡洛方法的出现得益于现代电子计算机的诞生，在1944年由Metropolis 和 Ulam提出于二战时美国原子弹研究的曼哈顿工程之中。蒙特卡洛这个名字是由Metropolis起的，借用了那个著名的赌场的名字，因为赌博总是和概率相关。