1
$\begingroup$

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:

enter image description here

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)
        
      answer_1 = num/denom
      return(answer_1)
    }

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)
            
          answer_1 = num/denom
          return(answer_1)
        }

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")

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.