问题描述
我需要使用以下内容创建随机游走算法 instructions
到目前为止，这是我所拥有的：
target = function(x){
return(ifelse(x<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 随机游走实验？
解决方法
创建频率分布
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()