###########################################################################################
#
# Random walk Metropolis-Hastings for the p-dimensional Gaussian target
# 
# The goal here is to study the behavior of single-move and block-move of the random 
# walk Metropolis-Hastings (RWMH) algorithm when the target distribution is a p-variate 
# Gaussian with zero mean and variance Sigma.  We study the behavior of the RWMH as 
# a function of the dimension p and the correlation rho.
#
# Hedibert F. Lopes
# 2026-05-19
#
###########################################################################################

# Log-density of N_p(0, Sigma) up to a constant
logtarget=function(theta){-0.5*as.numeric(t(theta)%*%SigmaInv%*%theta)}


p     = 10
rho   = 0.95

Sigma = matrix(rho, nrow = p, ncol = p)
diag(Sigma) = 1
SigmaInv = solve(Sigma)


# MCMC settings
burnin = 1000
lag    = 1
M      = 20000
sdev   = 1.0
niter  = burnin + lag*M
idx    = seq(burnin+1,niter,by=lag)


## Single-move 
theta       = rep(0, p)    
thetas.draw = NULL
accept_s    = 0
for (iter in 1:niter) {
  # update each component in turn
  for (j in 1:p) {               
    theta.star    = theta
    theta.star[j] = rnorm(1, theta[j], sdev)
    log_alpha     = logtarget(theta.star) - logtarget(theta)
    if (log(runif(1)) < log_alpha) {
      theta    = theta.star
      accept_s = accept_s + 1
    }
  }
  thetas.draw = rbind(thetas.draw, theta)
}
thetas.draw = thetas.draw[idx, ]

## Block-move 
theta        = rep(0, p)   
thetas.draw1 = NULL
accept_b     = 0
for (iter in 1:niter) {
	 # propose full vector
  theta.star = rnorm(p, theta, sdev)          
  log_alpha  = logtarget(theta.star) - logtarget(theta)
  if (log(runif(1)) < log_alpha) {
    theta    = theta.star
    accept_b = accept_b + 1
  }
  thetas.draw1 = rbind(thetas.draw1, theta)
}
thetas.draw1 = thetas.draw1[idx, ]

H     = 100
acfs  = matrix(0,H,p)
acfs1 = matrix(0,H,p)
for (i in 1:p){
  acfs[,i]=acf(thetas.draw[,i],plot=FALSE,lag.max=H)$acf[2:(H+1)]
  acfs1[,i]=acf(thetas.draw1[,i],plot=FALSE,lag.max=H)$acf[2:(H+1)]
}

par(mfrow=c(2,3))
ts.plot(thetas.draw,col=1:p,ylim=range(-5,5))
title(paste("(Burn-in,skip,M,sdev)=(",burnin,",",lag,",",M,",",sdev,")",sep=""))
ts.plot(acfs,ylim=c(0,1))
title(paste("Single-move acceptance rate:",round(accept_s/(niter*p),3),sep=""))
plot(density(thetas.draw[, 1]),main="",xlab="Components",xlim=c(-5,5))
for (i in 2:p)
  lines(density(thetas.draw[,i]),col=i)
curve(dnorm(x, 0, sqrt(Sigma[1, 1])), add = TRUE, lwd = 4)
title(paste("(p,rho)=(",p,",",rho,")",sep=""))

ts.plot(thetas.draw1,col=1:p,ylim=range(-5,5))
ts.plot(acfs1,ylim=c(0,1))
title(paste("Block-move acceptance rate:", round(accept_b/niter,3),sep=""))
plot(density(thetas.draw1[, 1]),main="",xlab="Components",xlim=c(-5,5))
for (i in 2:p)
  lines(density(thetas.draw1[,i]),col=i)
curve(dnorm(x, 0, sqrt(Sigma[1, 1])), add = TRUE, lwd = 4)



###########################################################################################
# Here we study the behavior of both RWMH algorithms (single-move vs block-move)
# for 9 combinations of the pair (p, rho)
###########################################################################################
ps     = c(2,10,50)
rhos   = c(0.5,0.75,0.95)

burnin = 10000
lag    = 1
M      = 10000
sdev   = 1.0
niter  = burnin + lag*M
idx    = seq(burnin+1,niter,by=lag)
theta.single = rep(0,niter)
theta.block  = rep(0,niter)
quants       = array(0,c(2,3,3,2))
for (ii in 1:3){
  for (jj in 1:3){
  	  p = ps[ii]
  	  rho = rhos[jj]
    Sigma = matrix(rho,p,p)
    diag(Sigma) = 1
    SigmaInv = solve(Sigma)
    theta       = rep(0,p)    
    for (iter in 1:niter){
      for (j in 1:p) {               
        theta.star    = theta
        theta.star[j] = rnorm(1, theta[j], sdev)
        log_alpha     = logtarget(theta.star) - logtarget(theta)
        if (log(runif(1)) < log_alpha) {
          theta    = theta.star
        }
      }
      theta.single[iter] = theta[1]
    }
    quants[1,ii,jj,] = quantile(theta.single,c(0.05,0.5))
    theta        = rep(0, p)   
    for (iter in 1:niter) {
      theta.star = rnorm(p, theta, sdev)          
      log_alpha  = logtarget(theta.star) - logtarget(theta)
      if (log(runif(1)) < log_alpha) {
        theta    = theta.star
      }
      theta.block[iter] = theta[1]
    }
    quants[2,ii,jj,] = quantile(theta.block,c(0.05,0.5))
  }
}

# Exact value is qnorm(0.5) = 0
tab1 = quants[1,,,2]
tab2 = quants[2,,,2]
rownames(tab1) = ps
colnames(tab1) = rhos
rownames(tab2) = ps
colnames(tab2) = rhos
round(cbind(tab1,tab2),3)

#        0.5   0.75   0.95    0.5   0.75  0.95
#  2  -0.034 -0.010 -0.033 -0.051 -0.022 0.029
# 10  -0.009  0.141  0.258  0.119  0.062 0.000
# 50  -0.006 -0.082 -0.446  0.000  0.000 0.000


# Exact value is qnorm(0.05) = -1.644854
tab1 = quants[1,,,1]
tab2 = quants[2,,,1]
rownames(tab1) = ps
colnames(tab1) = rhos
rownames(tab2) = ps
colnames(tab2) = rhos
round(cbind(tab1+1.644854,tab2+1.644854),3)

#        0.5   0.75   0.95    0.5   0.75  0.95
#  2  -0.027 -0.009 -0.018 -0.054 -0.041 0.027
# 10   0.002  0.014  0.314  0.200  0.363 1.645
# 50  -0.114 -0.155 -0.124  1.645  1.645 1.645