Gibbs Sampling in R

library(ggplot2)

Customize the 2D Gaussian Target

means <- c(x=1, y=1)
vars <- c(x=1, y=1)
corr <- -0.6
coordinates <- c("x", "y") # names of the coordinates, should match the names in the means and vars vectors.

Define the Marginal Distributions

Since the marginal distributions of x and y are both normal, a single function does the job.

# to sample the target p(x|y=1), use sample_marginal("x", 1)
sample_marginal <- function(sampled_coor = "x", other_value = 1){
  other_coor <- coordinates[coordinates != sampled_coor]  # "y" if sampled_var is "x"
  marginal_mean <- means[sampled_coor] + 
    corr * 
    sqrt(vars[sampled_coor]/vars[other_coor]) *
    (other_value - means[other_coor])
  
  marginal_var <- vars[sampled_coor] * (1-corr^2)
  
  return(rnorm(n=1, mean=marginal_mean, sd=sqrt(marginal_var)))
}

Define Each Iteration of Gibbs Sampling

# systematic scan Gibbs, simulates transition from prev to new 
sys_gibbs_step <- function(prev, order = coordinates){
  new <- prev
  for (coordinate in order){
    remaining_coord <- order[order != coordinate]  # the coord. not being updated
    new[coordinate] <- sample_marginal(sampled_var = coordinate,
                                      other_value = prev[remaining_coord])
  }
  return(new)
}

# random scan Gibbs
random_gibbs_step <- function(prev){
  new <- prev
  coordinate <- sample(coordinates, size=1) # sample the coordinate to be updated
  remaining_coor <- coordinates[coordinates != coordinate]
  new[coordinate] <- sample_marginal(sampled_coor = coordinate,
                                      other_value = prev[remaining_coor])
  return(new)
}

Running Gibbs Sampling

# initialize
n <- 50000  # takes a while... use a smaller sample size if your computer dies.
samples <- data.frame(x = numeric(length = n),
                      y = numeric(length = n))  # each realization is 2d.
previous_state <- c(x = xmean, y = ymean)

for (k in 1:n){
  previous_state <- random_gibbs_step(previous_state)  
  # previous_state <- sys_gibbs_step(previous_state)
  samples[k,] <- previous_state
}

Plotting the Samples

Beware that this first plot looks like crap if your sample size is small.

ggplot(data = samples, aes(x=x, y=y)) +
  geom_bin2d(bins = sqrt(n)) +
  scale_fill_continuous(type = "viridis") +
  theme_minimal()

So use this instead.

ggplot(data = samples, aes(x=x, y=y)) +
  geom_point(alpha = 5/sqrt(n)) +
  theme_minimal()