BMSI系列的第五篇,介绍贝叶斯计算中一个十分重要的算法——Metropolis Hastings

手敲$\LaTeX$难免出现纰漏,有任何疑似错误或者不清楚的地方请直接在下方评论区留言(评论区只开在【Series】下),谢谢各位读者。

M-H 算法

在高维空间中,上一节提到的简单采样方法往往不再奏效。此时我们使用另一个采样框架——MCMC。其核心思想是构造一个马尔可夫链Markov Cahin,其稳态分布为目标分布。然后通过在马尔可夫链上转移就完成了采样过程。本节提到的M-H算法 Metropolis-Hasting algorithm和下一节将提到的Gibbs sampling均属于这一框架。

下面我们记$p(\theta\lvert y)$为目标分布,$\theta^{(t)}$为$t$时刻我们采的样本。 M-H算法的框架如下

  1. 从提议函数(proposal)$g(\theta\lvert\theta^{(t)})$中获取一个提议$\theta^\star$(相当于前面的采样分布)

  2. 以概率$\min\{1,r\}$接受这个提议($\theta^{(t+1)}=\theta^\star$),否则拒绝这个提议($\theta^{(t+1)}=\theta^{(t)}$)。其中

$$r=r(\theta^{(t)},\theta^\star)=\frac{p(\theta^\star\lvert y)/g(\theta^\star\lvert\theta^{(t)})}{p(\theta^{(t+1)}\lvert y)/g(\theta^{(t)}\lvert\theta^\star)}=\frac{p(\theta^\star\lvert y)}{p(\theta^{(t+1)}\lvert y)}\cdot\frac{g(\theta^{(t)}\lvert\theta^\star)}{g(\theta^\star\lvert\theta^{(t)})}$$

其证明我们放在下一节介绍。我们这里先考察两个特殊的提议函数——独立的提议函数(即满足$g(\theta\lvert\theta^{(t)}=g(\theta)$)和对称的提议函数(即满足$g(\theta\lvert\theta^{(t)}=g(\theta^{(t)}\lvert\theta$)

独立的提议函数

此时,$r$的表达式中第二项简化为$g(\theta^{(t)})/g(\theta^\star)$

例子

$Y\sim\mathcal{N}(\theta,1)$,而参数有先验$\theta\sim t_1(0,1)$,此时参数后验分布为

$$p(\theta\lvert y)\propto p(y\lvert\theta)p(\theta)\propto \frac{\exp\{-(y-\theta)^2/2\}}{1+\theta^2}$$

我们选择提议函数为$\mathcal{N}(y,1)$,不难化简得此时

$$r=\frac{1+(\theta^{(t)})^2}{1+(\theta^\star)^2}$$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
MCMC <- function(n, burnout=0, start=0, seed=42) {
set.seed(seed)
N = n + burnout
samples = rep(NA, N)
current = start
for (i in 1:N) {
# set your proposal
proposal = rnorm(1, 1, 1) # y = 1
# set your ratio
logr = log(1 + current^2) - log(1 + proposal^2)
if (log(runif(1)) < logr) current = proposal
samples[i] = current
}
return(samples[(burnout+1):N])
}
N.factor <- integrate(function(x){exp(-(1-x)^2/2)/(1+x^2)}, -Inf, Inf)$value
posterior <- Vectorize(function(x){exp(-(1-x)^2/2)/(1+x^2) / N.factor})
hist(MCMC(100000), prob=T, breaks=30)
curve(posterior, -2, 4, add=T)
par(mfrow=c(2,1))
plot(x=1:100, y=MCMC(100), pch=20, main="x_0=0")
plot(x=1:100, y=MCMC(100, start=1000), pch=20, main="x_0=1000")
Independence_MH_1.png Independence_MH_1.png

重尾

回忆之前谈到的,拒绝采样法和重要性采样法都要求采样函数更加重尾,这里也不例外。否则,一旦$\theta$落入尾部,有$p(\theta^{(t)}\lvert y)\gg g(\theta^{(t)})$,对于绝大多数提议$\theta^\star$,有$p(\theta^\star\lvert y)\approx g(\theta^\star)$,此时

$$r=\frac{g(\theta^{(t)})}{p(\theta^{(t)}\lvert y)}\frac{p(\theta^\star\lvert y)}{g(\theta^\star)}\approx 0$$

也就是说绝大多数采样都会被拒绝,我们的转移陷入了停滞。让我们举一个重尾的例子:后验分布为$t_1$,而proposal使用$\mathcal{N}(0,1)$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
MCMC <- function(n, burnout=0, start=0, seed=42) {
set.seed(seed)
N = n + burnout
samples = rep(NA, N)
current = start
for (i in 1:N) {
# set your proposal
proposal = rnorm(1)
# set your ratio
logr = dnorm(current, log=T) - dnorm(proposal, log=T) + dt(proposal, 1, log=T) - dt(current, 1, log=T)
if (log(runif(1)) < logr) current = proposal
samples[i] = current
}
return(samples[(burnout+1):N])
}
par(mfrow=c(3,1))
plot(x=1:1000, y=MCMC(1000), pch=20, main="x_0=0")
plot(x=1:1000, y=MCMC(1000, start=5), pch=20, main="x_0=5")
plot(x=1:1000, y=MCMC(1000, start=4, burnout=10000), pch=20, main="x_0=5, burnout=10000")

可以看到,一旦采样到$(-3,3)$以外时,很容易就出现几十次的连续拒绝,这样的采样效果是很差的。而当初值极端时,甚至出现样本全部停滞的情况。一种对策是加入burnout的环节,使用一定的样本量进行“热身”,最后从“热身”后的起点开始抽样。但这也无法避免后续采样过程中出现极端样本后带来的影响,因此更明智的做法自然是仔细考察目标函数的性质,选用更加重尾的分布作为proposal。

对称的提议函数

此时,$r=\frac{q(\theta^\star\lvert y)}{q(\theta^{(t)}\lvert y)}$,此时这通常也被称为随机游走的Metropolis算法。

根据$r$的形式,不难看出,这个算法的特性是,它会接受所有的“更好的提议”,而有一定概率接受“更差的提议”

例子

仍使用上面的例子,只是proposal不再使用$\mathcal{N}(y,1)$这一独立函数,而是使用$\mathcal{N}(\theta^{(t)}, v^2)$这一对称函数。在这里,我们先设置$v^2=1$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
MCMC <- function(n, burnout=0, start=0, seed=42) {
set.seed(seed)
N = n + burnout
samples = rep(NA, N)
current = start
for (i in 1:N) {
# set your proposal
proposal = rnorm(1, current, 1)
# set your ratio
logr = -(1 - proposal)^2 / 2 + (1 - current)^2 / 2 + log(1 + current^2) - log(1 + proposal^2)
if (log(runif(1)) < logr) current = proposal
samples[i] = current
}
return(samples[(burnout+1):N])
}
N.factor <- integrate(function(x){exp(-(1-x)^2/2)/(1+x^2)}, -Inf, Inf)$value
posterior <- Vectorize(function(x){exp(-(1-x)^2/2)/(1+x^2) / N.factor})
hist(MCMC(100000), prob=T, breaks=30)
curve(posterior, -2, 4, add=T)
par(mfrow=c(2,1))
plot(x=1:100, y=MCMC(100), pch=20, main="x_0=0")
plot(x=1:100, y=MCMC(100, start=1000), pch=20, main="x_0=1000")
Randomwalk_MH_1.png Randomwalk_MH_2.png

调参

可以看到该方法也正确地完成了采样的工作。但当起始点偏离较多时,可以看到其采样的结果并非我们想要的——采样结果从1000缓慢的下降。可以想象,当我们设置一个较大的burnout后,采样工作可以步入正轨。但这仍是不让人满意的。为此,我们重新审视一下我们的采样过程。

我们每次的采样都是在前一步的基础上,在其周围选择下一个采样。而$v^2$控制的就是跳转的范围。当$v^2\approx 0$时,不难想象,每次的转移都只有极微小的变化,而此时的接受率$r\approx 1$,因此我们总是接受这些样本——于是我们得到了一大堆集中在某一个小区域的样本,这显然不是一个好的采样。而且当起始值较差时,就会出现样本极其缓慢地往正确的区域移动的情形(即上面出现的)。

当$v^2$很大时又会怎样呢?显然此时$\theta^\ast$很容易就跑到很极端的值,也就意味着$q(\theta^\ast\lvert y)\approx 0$,从而$r\approx 0$。此时绝大多数提议都被拒绝,采样效率也是极低的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
RW.MCMC <- function(n, sigma, burnout=0, start=0, seed=42) {
set.seed(seed)
N = n + burnout
samples = rep(NA, N)
current = start
for (i in 1:N) {
# set your proposal
proposal = rnorm(1, current, sigma)
# set your ratio
logr = -(1 - proposal)^2 / 2 + (1 - current)^2 / 2 + log(1 + current^2) - log(1 + proposal^2)
if (log(runif(1)) < logr) current = proposal
samples[i] = current
}
return(samples[(burnout+1):N])
}
par(mfrow=c(3,1))
plot(x=1:1000, y=RW.MCMC(1000, sigma=0.01), pch=20, main="sigma=0.01")
plot(x=1:1000, y=RW.MCMC(1000, sigma=1), pch=20, main="sigma=1")
plot(x=1:1000, y=RW.MCMC(1000, sigma=100), pch=20, main="sigma=100")
plot(x=1:1000, y=RW.MCMC(1000, start = 100, sigma=0.01), pch=20, main="sigma=0.01 start=100")
plot(x=1:1000, y=RW.MCMC(1000, start = 100, sigma=1), pch=20, main="sigma=1 start=100")
plot(x=1:1000, y=RW.MCMC(1000, start = 100, sigma=100), pch=20, main="sigma=100 start=100")
plot(x=1:1000, y=RW.MCMC(1000, start = 100, burnout=300, sigma=0.01), pch=20, main="sigma=0.01 start=100 burnout=300")
plot(x=1:1000, y=RW.MCMC(1000, start = 100, burnout=300, sigma=1), pch=20, main="sigma=0.01 start=100 burnout=300")
plot(x=1:1000, y=RW.MCMC(1000, start = 100, burnout=300, sigma=100), pch=20, main="sigma=0.01 start=100 burnout=300")
tuning_1.png tuning_2.png tuning_3.png

因此,应当存在一个合适的$v^2$来平衡大的跨度和高的接受率。然而手动尝试不同的$v^2$总不是一个好的解决方案。好在已有的一些研究表明,对于正态的proposal,最优的方差为$2.4^2\mathrm{Var}(\theta\lvert y) / d$,其中$d$是参数$\theta$的维度。而这个最优的方差会使得接受率约为$40\%$($d=1$),且随着参数维度的增加逐渐下降到$20\%$($d\rightarrow\infty$)。

那我们怎么获得$\mathrm{Var}(\theta\lvert y)$呢?能计算自然最好,如果难以计算,我们自然是使用抽样的方式来得到这个值的一个近似——只需要采出后验分布的样本然后计算方差即可。但这就牵扯到先有鸡还是先有蛋的问题了——我们求出这个值就是为了从后验分布采样而服务的,现在又要使用这个采样来求出这个值。因此,我们需要迭代式地进行这一过程。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
RW.Adapt.MCMC <- function(n, B, M, start=0, seed=42) {
set.seed(seed)
N = n + B * M
samples = rep(NA, N)
current = start
sigma = 1
for (b in 1:B) {
for (m in 1:M) {
i = (b - 1) * M + m
# set your proposal
proposal = rnorm(1, current, sigma)
# set your ratio
logr = -(1 - proposal)^2 / 2 + (1 - current)^2 / 2 + log(1 + current^2) - log(1 + proposal^2)
if (log(runif(1)) < logr) current = proposal
samples[i] = current
}
sigma = 2.4 * sd(samples[((b-1)*M):(b*M)])
}
for (i in (B*M+1):N) {
# set your proposal
proposal = rnorm(1, current, sigma)
# set your ratio
logr = -(1 - proposal)^2 / 2 + (1 - current)^2 / 2 + log(1 + current^2) - log(1 + proposal^2)
if (log(runif(1)) < logr) current = proposal
samples[i] = current
}
return(samples[(B*M+1):N])
}
N.factor <- integrate(function(x){exp(-(1-x)^2/2)/(1+x^2)}, -Inf, Inf)$value
posterior <- Vectorize(function(x){exp(-(1-x)^2/2)/(1+x^2) / N.factor})
hist(RW.Adapt.MCMC(10000, 10, 100, 100), prob=T, breaks=30)
curve(posterior, -2, 4, add=T)
par(mfrow=c(2,1))
plot(x=1:1000, y=RW.Adapt.MCMC(1000, 10, 100, 0), pch=20)
plot(x=1:1000, y=RW.Adapt.MCMC(1000, 10, 100, 100), pch=20)
length(unique(RW.Adapt.MCMC(1000, 10, 100))) / 1000 # acceotance ratio 0.452
tuning_4.png tuning_5.png