参考书目指出,如果 q 是对称分布,则比率 q(x|y)/q(y|x) 变为 1,该算法称为 Metropolis。那是对的吗?
是的,这是正确的。Metropolis 算法是 MH 算法的一个特例。
“随机漫步”大都会(-黑斯廷斯)怎么样?它与其他两个有何不同?
在随机游走中,提案分布在每一步之后重新以链最后生成的值为中心。通常,在随机游走中,提议分布是高斯分布,在这种情况下,这种随机游走满足对称性要求,算法是 Metropolis。我想您可以执行具有不对称分布的“伪”随机游走,这会导致提案在偏斜的相反方向上漂移(左偏分布会偏向右侧的提案)。我不知道你为什么要这样做,但你可以而且这将是一个大都会黑斯廷斯算法(即需要额外的比率项)。
它与其他两个有何不同?
在非随机游走算法中,提议分布是固定的。在随机游走变体中,提案分布的中心在每次迭代时都会发生变化。
如果提案分布是泊松分布怎么办?
然后你需要使用 MH 而不仅仅是 Metropolis。大概这将是对离散分布进行采样,否则您将不想使用离散函数来生成您的建议。
在任何情况下,如果采样分布被截断或者您事先知道它的偏斜,您可能想要使用非对称采样分布,因此需要使用 Metropolis-hastings。
有人可以给我一个简单的代码(C、python、R、伪代码或任何你喜欢的)示例吗?
这里是大都会:
Metropolis <- function(F_sample # distribution we want to sample
, F_prop # proposal distribution
, I=1e5 # iterations
){
y = rep(NA,T)
y[1] = 0 # starting location for random walk
accepted = c(1)
for(t in 2:I) {
#y.prop <- rnorm(1, y[t-1], sqrt(sigma2) ) # random walk proposal
y.prop <- F_prop(y[t-1]) # implementation assumes a random walk.
# discard this input for a fixed proposal distribution
# We work with the log-likelihoods for numeric stability.
logR = sum(log(F_sample(y.prop))) -
sum(log(F_sample(y[t-1])))
R = exp(logR)
u <- runif(1) ## uniform variable to determine acceptance
if(u < R){ ## accept the new value
y[t] = y.prop
accepted = c(accepted,1)
}
else{
y[t] = y[t-1] ## reject the new value
accepted = c(accepted,0)
}
}
return(list(y, accepted))
}
让我们尝试使用它来采样双峰分布。首先,让我们看看如果我们对 propsal 使用随机游走会发生什么:
set.seed(100)
test = function(x){dnorm(x,-5,1)+dnorm(x,7,3)}
# random walk
response1 <- Metropolis(F_sample = test
,F_prop = function(x){rnorm(1, x, sqrt(0.5) )}
,I=1e5
)
y_trace1 = response1[[1]]; accpt_1 = response1[[2]]
mean(accpt_1) # acceptance rate without considering burn-in
# 0.85585 not bad
# looks about how we'd expect
plot(density(y_trace1))
abline(v=-5);abline(v=7) # Highlight the approximate modes of the true distribution
现在让我们尝试使用固定的提案分布进行抽样,看看会发生什么:
response2 <- Metropolis(F_sample = test
,F_prop = function(x){rnorm(1, -5, sqrt(0.5) )}
,I=1e5
)
y_trace2 = response2[[1]]; accpt_2 = response2[[2]]
mean(accpt_2) # .871, not bad
起初这看起来不错,但如果我们看一下后密度......
plot(density(y_trace2))
我们会看到它完全停留在局部最大值。这并不完全令人惊讶,因为我们实际上将我们的提案分布集中在那里。如果我们将其集中在另一种模式上,也会发生同样的事情:
response2b <- Metropolis(F_sample = test
,F_prop = function(x){rnorm(1, 7, sqrt(10) )}
,I=1e5
)
plot(density(response2b[[1]]))
我们可以尝试在两种模式之间放弃我们的建议,但我们需要将方差设置得非常高,以便有机会探索其中的任何一种
response3 <- Metropolis(F_sample = test
,F_prop = function(x){rnorm(1, -2, sqrt(10) )}
,I=1e5
)
y_trace3 = response3[[1]]; accpt_3 = response3[[2]]
mean(accpt_3) # .3958!
请注意,我们的提案分布中心的选择如何对我们的采样器的接受率产生重大影响。
plot(density(y_trace3))
plot(y_trace3) # we really need to set the variance pretty high to catch
# the mode at +7. We're still just barely exploring it
我们仍然陷入两种模式中的较近者。让我们尝试在两种模式之间直接删除它。
response4 <- Metropolis(F_sample = test
,F_prop = function(x){rnorm(1, 1, sqrt(10) )}
,I=1e5
)
y_trace4 = response4[[1]]; accpt_4 = response4[[2]]
plot(density(y_trace1))
lines(density(y_trace4), col='red')
最后,我们越来越接近我们正在寻找的东西。从理论上讲,如果我们让采样器运行足够长的时间,我们可以从任何这些提议分布中得到一个有代表性的样本,但是随机游走非常快地产生了一个可用的样本,我们必须利用我们对后验如何假设的知识寻求调整固定的采样分布以产生可用的结果(说实话,我们还没有完全做到这一点y_trace4
)。
稍后我将尝试以大都会黑斯廷斯的示例进行更新。您应该能够相当容易地看到如何修改上述代码以生成 Metropolis Hastings 算法(提示:您只需将补充比率添加到logR
计算中)。