MCMC:Method 2 : MH Random walk with normal increments
require(mvtnorm)
mu <- c(1,2)
cov_mat <- matrix(c(1, 0.9, 0.9, 1), 2, 2, TRUE)
set.seed(1)
sample_chol <- rmvnorm(2000, mu, cov_mat)
x1<-sample_chol[,1]
y1<-sample_chol[,2]
sample_val <- function(x) {
delta_1 <- 0.6
delta_2 <- 0.4
return(x + c(rnorm(1, 0, sqrt(delta_1)), rnorm(1, 0, sqrt(delta_2))))
}
results_2 <- matrix(NA, nrow = N, ncol = 2)
results_2[1,] <- x_start
i <- 2
for( i in 2:N ) {
new_val <- sample_val(results_2[i - 1, ])
alpha_x_y <- alpha(results_2[i - 1], new_val)
if(runif(1) < alpha_x_y){
results_2[i, ] <- new_val
}else{
results_2[i,] <- results_2[i-1, ]
}
}
results_2 <- results_2[round(0.1*N):N,]
x<-results_2[,1]
y<-results_2[,2]
plot(
x,
y,
type=“p”,#线型
main=“cholesky”,W
sub=“MH???Random Walk”,
xlab=“x”,
ylab=“y”,
ylim=c(-6,6),
xlim=c(-6,6),
cex=0.3,#符号的大小
col="#cccccc",#必须要写全,六位的颜色标记;
pch=19,#设置标点的样式,19号表示实心;
asp=0.5)
#y/x的比例,y轴数值长度与x轴数值长度的比值
legend(“topleft”,inset=.02,c(“cholesky”,“MH-Random Walk”),pch=c(20,20),col=c(“red”,“grey”))
points(x1 , y1 ,pch=19, cex = 0.5, col= “red”)