一个简单的Metropolis算法的示例

Question

我需要用下面的instructions创建随机游走算法

到目前为止，我得到的是：

target = function(x){
  return(ifelse(x<0,0,exp(-x)))
}
x = rep(0,1000)
x[1] = 4     #initialize; I've set arbitrarily set this to 3
for(i in 2:10000){
  current_x = x[4]
  proposed_x = current_x + rnorm(1,mean=0,sd=1)
  A = target(proposed_x)/target(current_x) 
  if(runif(1)<A){
    x[i] = proposed_x       # accept move with probability min(1,A)
  } else {
    x[i] = current_x        # otherwise "reject" move, and stay where we are
  }
}
plot(x,main="values of x visited by the MH algorithm")
hist(x,xlim=c(0,10),probability = TRUE, main="Histogram of values of x visited by MH algorithm")
xx = seq(0,10,length=100)
lines(xx,target(xx),col="red")

我的代码的结果是：results 1& results 2。

如何重新创建MCMC随机游走实验？

Answer 1

创建频率分布

set.seed(8)
num_days <- 5e4
positions <- rep(0, num_days)
current <- 4
for (i in 1:num_days) {
# record current position
positions[i] <- current
# flip coin to generate proposal
proposal <- current + sample(c(-1, 1), size = 1)
# now make sure he loops around from 7 back to 1
if (proposal < 1) proposal <- 7
if (proposal > 7) proposal <- 1
# move?
prob_accept_the_proposal <- proposal/current
current <- ifelse(runif(1) < prob_accept_the_proposal, proposal, current)
}
tibble(theta = positions) %>%
ggplot(aes(x = theta)) +
geom_bar(fill = "steelblue") +
scale_x_continuous(expression(theta), breaks = 1:7) +
scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
theme_cowplot()

创建随机游走算法

tibble(t= 1:5e4,theta = positions) %>%
slice(1:500) %>%
ggplot(aes(x = theta, y = t)) +
geom_path(size = 1/4, color = "steelblue") +
geom_point(size = 1/2, alpha = 1/2, color = "steelblue") +
scale_x_continuous(expression(theta), breaks = 1:7) +
scale_y_log10("Time Step", breaks = c(1, 2, 5, 20, 100, 500)) + theme_cowplot()

创建参数直方图

tibble(x = 1:7, y = 1:7) %>%
ggplot(aes(x = x, y = y)) +
geom_col(width = .2, fill = "steelblue") +
scale_x_continuous(expression(theta), breaks = 1:7) +
scale_y_continuous(expression(p(theta)), expand = expansion(mult = c(0, 0.05))) +
theme_cowplot()