# Understanding the Metropolis-Hastings Algorithm

Can someone please tell me if I have correctly written the code for MCMC Sampling (Metropolis Hastings Algorithm)?

Given a multivariate joint probability distribution function, I am interested in learning how to randomly sample a conditional probability distribution function.

Suppose I have a 4 Dimensional multivariate Normal Distribution:

Suppose this multivariate Normal Distribution P(X, Y, Z, W) distribution has a:

Mean vector (4 x 1): 5.0060022 3.4280049 1.4620007 0.2459998

Variance-Covariance Matrix (4 x 4): 0.15065114 0.13080115 0.02084463 0.01309107 0.13080115 0.17604529 0.01603245 0.01221458 0.02084463 0.01603245 0.02808260 0.00601568 0.01309107 0.01221458 0.00601568 0.01042365

Using the R programming language, I created a function that corresponds to this 4 Dimensional Multivariate Normal Distribution:

#define constants needed for the multivariate normal
sigma1 <- c(0.15065114 , 0.13080115 ,  0.02084463 , 0.01309107 , 0.13080115 , 0.17604529 ,  0.01603245 , 0.01221458 , 0.02084463 , 0.01603245  , 0.02808260 , 0.00601568 , 0.01309107 , 0.01221458 ,  0.00601568 , 0.01042365)

sigma <- matrix(sigma1, nrow=4, ncol= 4, byrow = TRUE)
sigma_inv <- solve(sigma)
sigma_det <- det(sigma)
denom = sqrt( (2*pi)^4 * sigma_det)

#actual multivariate function is defined below ("target")
target <- function(x,y,z,w)

{
x_one = x - 5.0060022
x_two = y - 3.4280049
x_three = z - 1.4620007
x_four = w - 0.2459998

x_t = c(x_one, x_two, x_three, x_four)
x_t_one <- matrix(x_t, nrow=4, ncol= 1, byrow = TRUE)
x_t_two =  matrix(x_t, nrow=1, ncol= 4, byrow = TRUE)

# In this part, as it's (x-mu)^T * SIGMA * (x-mu)

#num = exp(-0.5 * t(x_t_t) %*% sigma1_inv %*% x_t_t)
num = exp(-0.5 * x_t_two  %*%  sigma_inv  %*%  x_t_one)

}


Question: Suppose I want to take random samples from this multivariate normal distribution P(X, Y, Z, W), conditional on P(X, Y | Z = 2 , W = 1.3)

I attempted to manually write a Monte Carlo Sampler (Metropolis-Hastings) to take random samples from P(X, Y | Z = 2 , W = 1.3). To do this, I first "fixed" the values of Z and W within the original multivariate normal distribution:

#fix the definitions of w and z as per P(X, Y | Z = 2 , W = 1.3)
target <- function(x,y)

{
x_one = x - 5.0060022
x_two = y - 3.4280049
x_three = 2 - 1.4620007
x_four = 1.3 - 0.2459998

x_t = c(x_one, x_two, x_three, x_four)
x_t_one <- matrix(x_t, nrow=4, ncol= 1, byrow = TRUE)
x_t_two =  matrix(x_t, nrow=1, ncol= 4, byrow = TRUE)

# In this part, as it's (x-mu)^T * SIGMA * (x-mu)

#num = exp(-0.5 * t(x_t_t) %*% sigma1_inv %*% x_t_t)
num = exp(-0.5 * x_t_two  %*%  sigma_inv  %*%  x_t_one)

}


Next, I run the Monte Carlo Sampler (Metropolis-Hastings) to randomly sample this Conditional Distribution:

library(mvtnorm)
x = rep(0,10000)
y = rep(0,10000)

x[1] = 3     #initialize; I've set arbitrarily set this to 3 and 1
y[1] =1

for(i in 2:10000){
current_x = x[i-1]
current_y = y[i-1]

new = rmvnorm(n=1, mean=c(current_x,current_y), sigma=diag(2), method="chol")   # generate bivariate random numbers
proposed_x = new[1]
proposed_y = new[2]

A = target(proposed_x,proposed_y)/target(current_x,current_y)
if(runif(1)<A){
x[i] = proposed_x       # accept move with probabily min(1,A)
y[i] = proposed_y

} else {
x[i] = current_x        # otherwise "reject" move, and stay where we are
y[i] = current_y

}
}


The final answer can be taken as the "mean" of both columns:

mean(mcmc_output$y) mean(mcmc_output$y)


My Question: Can someone please tell me if what I have done is correct? I know that the Multivariate Normal Distribution has "attractive theoretical properties" that allow you to perform the above task without using MCMC Sampling - but I am trying to learn how to perform MCMC Sampling and wanted to start with a basic example by using the Normal Distribution. Can someone please tell me if what I have done is correct? Is there a way to find out if the MCMC algorithm has converged?

Thanks!

Optional: Visualization of the Metropolis-Hastings Algorithm

library(ggplot2)

ggplot(mcmc_output, aes(x = x, y =
y)) +
geom_density_2d_filled() +
ggtitle("Contour Plots of the MCMC Estimates")