贝叶斯集锦（3）：从MC、MC到MCMC
2013-07-31 23:03:39
#####一份草稿
贝叶斯计算基础
一、从MC、MC到MCMC
斯坦福统计学教授Persi Diaconis是一位传奇式的人物。Diaconis14岁就成了一名魔术师,为了看懂数学家Feller的概率论著作，24岁时进入大学读书。他向《科学美国人》投稿介绍他的洗牌方法，在《科学美国人》上常年开设数学游戏专栏的著名数学科普作家马丁•加德纳给他写了推荐信去哈佛大学，当时哈佛的统计学家Mosteller 正在研究魔术，于是Diaconis成了Mosteller的学生。（对他这段传奇经历有兴趣的读者可以看一看统计学史话《女士品茶》）。 下面要讲的这个故事，是Diaconis 在他的文章The Markov Chain Monte Carlo Revolution中给出的破译犯人密码的例子。 一天，一位研究犯罪心理学的心理医生来到斯坦福拜访Diaconis。他带来了一个囚犯所写的密码信息。他希望Diaconis帮助他把这个密码中的信息找出来。 这个密码里的每个符号应该对应着某个字母，但是如何把这些字母准确地找出来呢？Diaconis和他的学生Marc采用了一种叫做MCMC（马尔科夫链蒙特卡洛）的方法解决了这个问题。这个MCMC方法就是这一节我们所要讨论的内容。
(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起的，借用了那个著名的赌场的名字，因为赌博总是和概率相关。
例：蒲丰投针 1777年，法国学者蒲丰提出的一种用随机试验来求圆周率的方法
set.seed(1234)
l <- 0.8 # 针的长度
n <- 1e+06 # 重复100000次.
u1 <- runif(n) #取随机数
x <- 1/2 * u1 # x是针中心到最近的线的距离
u2 <- runif(n)
y <- l/2 * sin(u2 * 2 * pi)
z <- as.numeric(x <= y) # 相交的充要条件是x<=y。
pi.e <- n * l/sum(z) # π的估计式
pi.e
## [1] 3.141
蒙特卡洛模拟的过程是随机的，但解决的不仅有随机的问题也可以有确定性的问题。蒙特卡洛可以视为一种思想的泛称，只要在解决问题分人过程中，利用大量的随机样本，然后对这些样本结果进行概率分析从而得到问题求解的方法，都可以称之为蒙特卡洛方法。 蒙特卡洛方法的基础是（利用计算机）生成指定分布的随机数的抽样。在统计上发展了一些方法来解决这个问题。一般来说，从均匀分布U（0,1）的样本较容易生成，通过线性同余发生器可以生成伪随机数序列，这些序列的各种统计指标和均匀分布 U(0,1) 的理论计算结果非常接近。这样的伪随机序列就有比较好的统计性质，可以被当成真实的随机数使用。常见的分布的随机数都可以基于均匀分布的样本生成。 对定积分的计算是蒙特卡洛方法的一种重要的应用，而很多统计问题，比如计算概率、矩（期望、方差等）都可以归结为定积分的计算比如说我们想计算θ的期望，期望的计算公式为：
E(θ)=∫θp(θ)dθ
怎么样通过抽样的方法来解决这个问题呢？
在θ的分布中抽取独立样本θ1,θ2,...,θn，n足够大
用样本均值θˉ=1n∑ni=1θi作为E(θ)
当n趋近无穷时，利用大数定律，样本均值会收敛到E(θ)（这称为仿真一致性）
从上面这个简单例子，蒙特卡洛方法的基本思路是：
（a）针对实际问题建立一个简单易行的概率统计模型，使问题所求的解为该模型的概率分布或者数字特征，比如：某个事件的概率或者是某个随机变量的期望值。
（b）对模型中的随机变量建立抽样方法，在计算机上进行模拟试验，得到足够的随机抽样，并对相关事件进行统计
（c）对试验结果进行分析，给出所求解的估计及其精度（方差）的估计
上面例子中所用的这种模拟方法是把有限区间上的积分看做这个区间上均匀分布随机变量的期望值，叫做逆分布函数法（inverse-CDF
method，也称作逆变换法）。这种方法简便易行，但对于无穷区间的积分是不适用的。所以还有一些其它的抽样方法，比如importance
sampling;rejection sampling;slice sampler等，不再一一详述。
对于复杂或者高维的分布，利用蒙特卡洛方法生成随机样本是比较困难的，因此需要发展一些更复杂的随机模拟技术。
（3）马尔科夫链（Markov Chain）和MCMC
在蒙特卡洛模拟中，我们在后验分布中抽取样本，当这些样本独立时，利用大数定律样本均值会收敛到期望值。如果得到的样本是不独立的，那么就要借助于马尔科夫链进行抽样。MCMC方法就是为了这个目的而诞生的。
马尔科夫链又称为马尔科夫过程是一种离散的随机过程，用俄罗斯数学家安德烈.马尔科夫的名字命名。随机过程可以看做一个随时间变化的随机变量序列，在物理上用来表示在某个空间中物体的位置随机变化的轨迹；在贝叶斯统计中，物体的轨迹可以看做蒙特卡洛算法的结果，而“空间”就是支持后验分布p(θ∣x)的样本空间。
马尔科夫链的定义如下， 设θ(t)是一个随机过程，如果它满足下面的性质： p(θ(t+h)=xt+h∣θ(s)=xs,s≤t)=p(θ(t+h)=xt+h∣θ(t)=xt,∨h>0
这个定义又称为马尔科夫性质，对一个马尔科夫链来说，未来状态只与当前t时刻有关，而与t时刻之前的历史状态无关（条件独立）。
在时刻t，一个遍历的马尔科夫链从状态i可以向链的任一状态（包括状态i自身）转移。从状态i到状态j的转移概率记作p(i,j)或者p(i\rightarrow j),p(j\mid i)。所有这些转移概率构成了状态转移矩阵，记为P。
马尔科夫链的一个很重要的性质是平稳分布。简单的说，主要统计性质不随时间而变的马尔科夫链就可以认为是平稳的。数学上有马氏链收敛定理，当步长n足够大时，一个非周期且任意状态联通的马氏链可以收敛到一个平稳分布π(x)。这个定理就是所有的MCMC方法的理论基础。
我们对一个马氏链进行状态转移。设从初始分布π0出发，假设到第n步这个链可以收敛到平稳分布π0，那个这个过程可以表示为： X0 π0(x),X1
π1(x),...,Xn πn(x)=π(x),Xn+1 π(x),Xn+2 π(x),...
在利用马氏链进行抽样时，在收敛之前的一段时间，比如上面的前n-1次迭代，各个状态的边际分布还不能认为是稳定分布，所以在进行估计的时候，应该把前面的这n-1个迭代值去掉。这个过程称为“burn-in”。
MCMC方法就是构造合适的马尔科夫链进行抽样而使用蒙特卡洛方法进行积分计算。 既然马尔科夫链可以收敛到平稳分布。我们可以建立一个以π为平稳分布的马尔科夫链，对这个链运行足够长时间之后，可以达到平稳状态。此时马尔科夫链的值就相当于在分布π(x)中抽取样本。利用马尔科夫链进行随机模拟的方法就是MCMC。这种方法最早是由 Metropolis研究物理学中粒子系统的平稳性质想到的，他得到了第一个MCMC的算法：Metropolis算法。由此，一系列的MCMC方法相继诞生。
二、M-H算法
MCMC方法最早由Metropolis（1954）给出，后来Metropolis的算法由Hastings改进，合称为M-H算法。M-H算法是MCMC的基础方法。由M-H算法演化出了许多新的抽样方法，包括目前在MCMC中最常用的Gibbs抽样也可以看做M-H算法的一个特例。
我们来介绍基本的M-H算法。马氏链的收敛性质主要由转移矩阵P 决定, 所以基于马氏链做抽样的关键问题是如何构造转移矩阵P,使得平稳分布恰好是我们要的分布π(x)
这里用到了马尔科夫链的另一个性质，如果具有转移矩阵P和分布π(x)的马氏链对所有的状态i,j满足下面的等式：
π(i)p(i,j)=π(j)p(j,i)
这个等式称为细致平衡方程。满足细致平衡方程的分布π(x)是平稳的。 所以我们希望抽样的马尔科夫链是平稳的，可以把细致平衡方程作为出发点。
但是一般情况下，任给的分布分布π(x)不一定是平稳的。那么为了让细致平衡方程成立，我们可以引入一个函数α(i,j) 满足≤α(i,j)≤1,使得
π(i)p(i,j)α(i,j)=π(j)p(j,i)α(j,i)
按照对称性，我们可以取：
α(i,j)=π(j)p(j,i)α(j,i)=π(j)p(i,j)
记
q(i,j)=p(i,j)α(i,ij)q(j,i)=p(j,i)α(j,i)
这样就得到了新的以q(i,j)为转移概率的马尔科夫链，且具有平稳分布π(x)
这里的α(i,j)叫做接受率，也就是在原来的马氏链上以概率α(i,j)接受从状态i到状态j的转移概率为p(i,j)的状态转移。 进一步的，如果接受率过小，会导致马氏链过多地拒绝状态转移，这样马氏链的收敛速度会很慢。因此，我们可以考虑放大α(i,j)，使 π(i)p(i,j)α(i,j)=π(j)p(j,i)α(j,i)中的两个方程中大的\alpha$取到1，小的同比例放大。这样拒绝率就可以表示为：
α(i,j)=min(1,π(j)p(j,i)/π(i)p(i,j))
这就是M-H算法。 我们的讨论都是针对离散状态的，对连续状态的马尔科夫链依然有相同的结果。在连续状态中表示状态转移概率的项用条件概率密度代替，称为转移核。
M-H算法的步骤：
(1)构造合适的提议分布（Proposal distribution) g(⋅∣Xt)在分布g中产生X0；
(2)迭代下面的步骤：
在g(⋅∣Xt)中生成新样本Y；
从均匀分布U(0,1)抽取随机数U；
如果U满足U≤f(Y)g(Xt∣Y)/f(Xt)g(Y∣Xt), 则令Xt+1=Y（转移到新状态），否则Xt+1=Xt（状态不变）。其中f是目标分布，也就是 我们需要进行抽样的后验分布；
增加t值，进行下一步迭代。
例：对自由度为v的t分布进行抽样， t分布的密度值由R函数dt来计算
N <- 3000 #抽样个数
v <- 4 #自由度
x <- c()
x0 <- 0 # 初始值
x[1] <- x0
k <- 0 #k表示拒绝转移的次数
u <- runif(N) #抽取均匀分布随机数
for (i in 2:N) {
y <- rnorm(1, x[i - 1], 2)
if (u[i] <= (dt(y, v)/dt(x[i - 1], v)))
x[i] <- y else {
x[i] <- x[i - 1]
k <- k + 1
}
}
k
## [1] 1361
回到我们前面的故事。因为在文本中前后字母有依赖的关系，所以Marc利用前后字母的依赖关系，来描述整个密码文本中出现这些字符的概率模型。他用《战争与和平》作为标准的文本，统计一个字母到另一个字母的一步转移概率。然后他利用了M-H算法，在假设所有对应关系出现的可能性相等的前提下(也就是无信息先验的情况)，首先随机给出了字符和字母的对应关系。利用前边得到的转移概率，计算这种对应关系出现的概率p1。随机抽取两个字符，互换它们的对应关系，此时对应概率p2。如果P2>P1，接受新的对应关系；否则，抛一枚以P2/P1的概率出现正面的硬币，如果出现正面，则接受新的对应关系，否则依然保持旧有的对应关系。这就是MH算法的运用。当算法收敛时， 就会得到真实的对应关系。事实上，当算法运行了2000多步的时候，就得到了有意义的信息，一个混合了英语和西班牙语的文本段落。
三、Gibbs抽样
Gibbs抽样是 Stuart Geman 和Donald Geman 这两兄弟于1984年提出来的,可以看做MH算法当α=1的一个特例，用于目标分布为多元分布的情况。
假设在多元分布中所有的一元条件分布（每个分量对所有其它分量的条件分布，这个分布也叫做满条件分布）都是可以确定的。记m维随机向量X=(X1,X2,...,Xm)′，X−i表示X中去掉分量Xi后剩余的m-1维向量。那么一元条件分布就是f(Xi∣X−i)。
Gibbs抽样就是在这m个条件分布中迭代产生样本。
算法步骤如下： (1)给出初值X(0)； (2)对t=1,…,T，进行迭代： (a)令x1=X1(t−1)； (b)依次更新每一个分量，即对i=1,…,m, (i)从f(Xi∣X−i(t−1))中产生抽样Xi(t)； (ii)更新xi(t)=Xi(t)； (c)令X(t)=(X1(T),X2(T),...,Xm(T))′(每个抽取的样本都被接受了)； (d)更新t。
在这个算法里，对每一个状态t,X(t)的分量是依次更新的。这个分量更新的过程是在一元分布f(Xi∣X−i)中进行的，所以抽样是比较容易的。
例：使用Gibbs抽样抽取二元正态分布N(μ1,μ2,σ12,σ22,ρ)的随机数 在二元正态分布的条件下，两个分量的一元条件分布依然是正态分布：
f(x1∣x2)∼N(μ1+ρσ1σ2(x2−μ2),(1−ρ2)σ12)f(x2∣x1)∼N(μ2+ρσ2σ1(x1−μ1),(1−ρ2)σ22)
采用上面所述的Gibbs抽样方法，对给定参数值的一个二元正态分布进行抽样。
n <- 5000 #抽样个数（链的长度）
burn.in <- 2500 #前2000个抽样按burn-in处理
X <- matrix(0, n, 2)
mu1 <- 1 #对参数赋值
mu2 <- -1
sigma1 <- 1
sigma2 <- 2
rho <- 0.5
s1.c <- sqrt(1 - rho^2) * sigma1
s2.c <- sqrt(1 - rho^2) * sigma2
X[1, ] <- c(mu1, mu2) #初始化
for (i in 2:n) {
x2 <- X[i - 1, 2]
m1.c <- mu1 + rho * (x2 - mu2) * sigma1/sigma2
X[i, 1] <- rnorm(1, m1.c, s1.c)
x1 <- X[i, 1]
m2.c <- mu2 + rho * (x1 - mu1) * sigma2/sigma1
X[i, 2] <- rnorm(1, m2.c, s2.c)
}
b <- burn.in + 1
x.mcmc <- X[b:n, ]
library(MASS)
MVN.kdensity <- kde2d(x.mcmc[, 1], x.mcmc[, 2], h = 5) #估计核密度