BMSI系列的第六篇,介绍贝叶斯计算中介绍用于多维分布中采样的Gibbs Sampling

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

Gibbs Sampling

Gibbs Sampling的思想是,将高维空间的采样转化为若干个一维空间的采样,而这是我们所熟悉的。具体的,对于,我们的采样过程为

当然,事实上每个也并不要求是一个标量/长度为1的向量。如果你觉得联合抽样更方便(比如是二元正态分布)也可以一次抽出。此时也叫做分块的Gibbs采样(block Gibbs sampling)

二元正态分布

我们首先用一个简单的例子来展示Gibbs Sampling。假设我们的目标分布为。根据概率论知识我们不难求得,于是乎

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
Gibbs.example1 <- function(n, rho, start=c(0,0), seed=42) {
set.seed(seed)
sigma = sqrt(1 - rho^2)
samples = matrix(NA, nrow=n, ncol=2)
samples[1,] = start
for (i in 2:n) {
samples[i, 1] = rnorm(1, rho*samples[i-1, 2], sigma)
samples[i, 2] = rnorm(1, rho*samples[i, 1], sigma)
}
return(samples)
}
samples = Gibbs.example1(1000, 0.8)
plot(x=samples[,1], y=samples[,2], pch=20)
colMeans(samples)
# [1] -0.04617112 -0.04515437
cov(samples)
# [,1] [,2]
# [1,] 1.0121685 0.8152691
# [2,] 0.8152691 1.0126826

Gibbs + M-H

可以看到,Gibbs抽样的核心就是从中进行抽样(表示除去后的所有)然而有时候从这个分布中直接抽样并不容易,于是我们可以嵌入M-H算法完成这一抽样。仍然使用上面的例子,但我们假装没办法完成这一抽样,于是使用Randomwalk MH进行采样

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
Gibbs.example1.2 <- function(n, rho, start=c(0,0), seed=42) {
set.seed(seed)
sigma = sqrt(1 - rho^2)
samples = matrix(NA, nrow=n, ncol=2)
samples[1,] = start
for (i in 2:n) {
samples[i, 1] = rnorm(1, rho*samples[i-1, 2], sigma)
proposal = rnorm(1, samples[i-1, 2], 2.4*v) # Notice that 2.4*v is the optimal proposal sd
logr = dnorm(proposal, rho*samples[i, 1], v, log=T)
- dnorm(samples[i-1, 2], rho*samples[i, 1], v, log=T)
samples[i, 2] = ifelse(log(runif(1))<logr, proposal, samples[i-1, 2])
}
return(samples)
}
samples = Gibbs.example1.2(1000, 0.8)
plot(x=samples[,1], y=samples[,2], pch=20)
colMeans(samples)
# [1] -0.04617112 -0.04515437
cov(samples)
# [,1] [,2]
# [1,] 1.0121685 0.8152691
# [2,] 0.8152691 1.0126826

参数扩张

参数扩张的思想是引入额外的参数,而目标分布则是扩张后的分布的一个边缘分布。我们通过在这个扩张后的空间内采样,也就得到了目标分布的采样。为方便计,我们假设目标的参数是一个标量,目标分布为。注意到

因此可以视为下面这个均匀分布的一个边缘分布

因此我们的采样过程为

我们使用一个简单的例子,假设,即,此时

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Gibbs.example2 <- function(n, start=1, seed=42) {
set.seed(seed)
u = rep(NA, n)
samples = rep(NA, n)
samples[1] = start
for (i in 2:n) {
u[i] = runif(1, 0, exp(-samples[i-1]))
samples[i] = runif(1, 0, -log(u[i]))
}
return(samples)
}
hist(Gibbs.example2(10000), prob=T, breaks=30)
curve(exp(-x), 0, 8, add=T)

更多的例子

Hierarchical normal example

模型假设为,其中。我们选取先验

可以推导出

根据上式既可以完成Gibbs采样。

我们还可以选取另一个先验

其中为柯西分布只取正半轴的分布(乘二)。得到的形式与之前几乎一样,除了

此时从抽样就有许多方法了,比如

  • 使用独立的M-H算法,从中采样出,然后以接受率接受这一转移
  • 使用随机游走的M-H算法,使用某个对称的proposal(如正态),然后以接受率接受这一转移

Hierarchical binomial example

模型假设为。我们选取先验

可以推导出

注意这里的的抽样彼此是独立的,因此可以一起抽取,算是 block Gibbs sampling里一个特殊情形。

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
39
40
log_alpha = function(theta, alpha, beta) {
if (alpha < 0) return(-Inf)
n = length(theta)
return((alpha - 1) * sum(log(theta))
- n * lbeta(alpha, beta) - 5 / 2 * (alpha + beta))
}
log_beta = function(theta, alpha, beta) {
if (beta<0) return(-Inf)
n = length(theta)
return((beta - 1) * sum(log(1 - theta))
- n * lbeta(alpha, beta) - 5 / 2 * (alpha + beta))
}
MCMC = function(n, dat, inits, tune) {
n_groups = nrow(dat)
alpha = inits$alpha
beta = inits$beta
samples = matrix(NA, nrow=n, ncol=2+n_groups)
for (i in 1:n) {
# Sample thetas
theta = with(dat, rbeta(length(y), alpha+y, beta+n-y))
# Sample alpha
alpha_prop = rnorm(1, alpha, tune$alpha)
logr = log_alpha(theta, alpha_prop, beta)-log_alpha(theta, alpha, beta)
alpha = ifelse(log(runif(1))<logr, alpha_prop, alpha)
# Sample beta
beta_prop = rnorm(1, beta, tune$beta)
logr = log_beta(theta, alpha, beta_prop)-log_beta(theta, alpha, beta)
beta = ifelse(log(runif(1))<logr, beta_prop, beta)
# Record parameter values
samples[i,] = c(alpha, beta, theta)
}
colnames(samples) <- c('alpha', 'beta', paste('theta', 1:71, sep='_'))
return(samples)
}
rats <- read.table("rats_data.txt",header=T)
inits = list(alpha=1, beta=1)
samples = MCMC(2000, dat=rats, inits=inits, tune=list(alpha=1,beta=1))

其中数据文件在层次化模型那一节出现过,为方便这里也放上链接数据