I am working with the R programming language. I have been trying to work on the following problem: Given a multivariate joint probability distribution function, I am interested in learning how to randomly sample (using MCMC) a conditional probability distribution function.
At the moment, I have found no "pre-made" way to do this, and have had to write the R code myself. The R code that I have written seems to work - but I am not sure if the R code is correct. I would be interested in finding a "pre-made" (e.g. some R package/library) which can do this. Below, I show the steps I have taken so far to solve this problem:
Part 1 - Context: 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):
[2,] 0.13080115 0.17604529 0.01603245 0.01221458
[3,] 0.02084463 0.01603245 0.02808260 0.00601568
[4,] 0.01309107 0.01221458 0.00601568 0.01042365 ```
Part 2 - Function Definition
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)
}
Part 3 - MCMC Code:
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 results of the MCMC can be seen by calling [x] and [y]. Thus, the final estimates would be:
mean(mcmc_output$x)
[1] 6.281715
mean(mcmc_output$y)
[1] 4.63817
(Optional) Part 4 - Visualize Results:
Now, I can plot the histograms for the results:
mcmc_output = data.frame(x,y)
par(mfrow=c(1,2))
plot(density(mcmc_output$y, main = "Density of Y"))
plot(density(mcmc_output$x, main = "Density of X"))
I can also look at the contour plots of these MCMC samples:
library(ggplot2)
ggplot(mcmc_output, aes(x = x, y =
y)) +
geom_density_2d_filled() +
ggtitle("Contour Plots of the MCMC Estimates")
Question: Is there a "easier" (standard) way to do this (MCMC Sampling of Custom Defined Functions) in R using a well known and pre-built package/library? The code I have written seems long and inefficient - and probably contains some mistakes as well. Does anyone know about, or can anyone recommend a package/library in R that does this?
Thanks!
References: