对于部分概率分布,无法采用积分的方法计算其期望(特别是在贝叶斯后验分布中,我们算不出来分母上的积分项)这时候可以采用MCMC
的方法来进行采样,基于采样样本对分布情况进行估计。
1. Metropolis-Hastings方法
我们的目标函数是,初始采样值为
。我们可以随机设定一个转移矩阵
,即
,然后抽取随机数,每次抽取按照
的概率接受或拒绝新抽取的数。若接受,
,否则,
。
2. 案例及采样逻辑
例如,我们的目标函数是混合正态分布:
我们首先定义转移矩阵,也就是可以任意指定的。我们定义为:
其中, 是一个以零为中心的对称分布
,这样就有:
而且:
这样,我们的接受概率就可以写成:
3. 实操
- 定义目标函数:
f <- function(x) {
0.3 * dnorm(x, 0, 1) + 0.7* dnorm(x, 4, 1)
}
- 任意指定
,即
,然后让
g <- function(x) rnorm(1, x, 4)
- 然后我们定义具体的抽样过程:
# 初始值;目标函数f;指定的任意函数g;抽样次数
rw.mh <- function(x0, f, g, n.steps = 1000, ...){
x <- x0
samples <- matrix(NA, n.steps, length(x))
for (i in seq_len(n.steps)){
# 猜想
x.new <- g(x)
# 计算接受概率
alpha <- min(1, f(x.new, ...) / f(x, ...))
# 判断是否接受
if (runif(1) < alpha){
x <- x.new
}
samples[i,] <- x
}
return(samples)
进行抽样:
x0 <- 1
n.steps <- 5000
samples <- rw.mh(x0, f, g, n.steps = n.steps)
我们最后比较一下,实际的概率密度函数,和抽取样本的概率密度函数的情况:
p1 <- ggplot(data.frame(x=runif(3000, -5, 15)), aes(x)) +
stat_function(fun = f) +
ylab('f(x)')
# 实际的
p2 <- ggplot(mapping = aes(x = samples[,1])) +
geom_density() +
ylab('f(x)')
Fig.1