Metropolis 采样

对于部分概率分布,无法采用积分的方法计算其期望(特别是在贝叶斯后验分布中,我们算不出来分母上的积分项)这时候可以采用MCMC的方法来进行采样,基于采样样本对分布情况进行估计。

1. Metropolis-Hastings方法

我们的目标函数是f(x),初始采样值为x_0。我们可以随机设定一个转移矩阵Q,即g\left(x^{*} \mid x_{t-1}\right),然后抽取随机数,每次抽取按照
\alpha=\min \left(1, \frac{f\left(x^{*}\right) g\left(x_{t-1} \mid x^{*}\right)}{f\left(x_{t-1)} g\left(x^{*} \mid x_{t-1}\right)\right)}\right)
的概率接受或拒绝新抽取的数x^*。若接受, x_t = x^*,否则,x_t=x_{t-1}

2. 案例及采样逻辑

例如,我们的目标函数是混合正态分布:
f(x) \sim 0.3 * N(0,1) + 0.7*N(4,1)

我们首先定义转移矩阵,也就是可以任意指定的g\left(x^{*} \mid x_{t-1}\right)。我们定义为:

x^{*}=x_{t-1}+\epsilon

其中,\epsilon 是一个以零为中心的对称分布q,这样就有:

g\left(x^{*} \mid x_{t-1}\right)=q(\epsilon)
而且:
g\left(x_{t-1} \mid x^{*}\right)=q(-\epsilon)=q(\epsilon)

这样,我们的接受概率就可以写成:

\alpha=\min \left(1, \frac{f\left(x^{*}\right)}{f\left(x_{t-1}\right)}\right)

3. 实操

  • 定义目标函数:
f <- function(x) {
    0.3 * dnorm(x, 0, 1) + 0.7* dnorm(x, 4, 1)
}
  • 任意指定g\left(x^{*} \mid x_{t-1}\right),即x^{*}=x_{t-1}+\epsilon,然后让\epsilon \sim N(0,4^2)
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
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。