R/update_sigma.ab_horseshoe_fc.R

In ianmtaylor1/amen2: Alternate MCMC Algorithm for AME Models

# Function to calculate the log density of the full conditional distribution
# of sigma.a or sigma.b (up to porportionality). Argument is called 'a' but it's
# the same function for b, also
# x - value of sigma.a (or sigma.b)
# a - vector of values of a (or b)
sigma.ab.fc.log <- function(x, a) {
  n <- length(a)
  return(-n * log(x) - sum(a^2) / (2 * x^2) - log(1 + x^2))
}

# Simulates from the full conditional distribution for sigma.ab
# Under the global horseshoe prior. Uses Metropolis-Hastings and depends
# On the current value of sigma.ab
update_sigma.ab_horseshoe_fc <- function(Y, Z, X, beta, sigma.squared, a, b, sigma.ab) {
  # Number of times to "jump"
  iter <- 25
  # Jump kernel - x is the current value to jump from
  jump <- function(x) { return( abs(rnorm(n=1, mean=x, sd=2)) ) }
  # Current values of sigma.a and sigma.b
  sigma.a <- sqrt(sigma.ab[1,1])
  sigma.b <- sqrt(sigma.ab[2,2])
  # Current values of the "density" for determining jump probabilities
  current.logf.a <- sigma.ab.fc.log(sigma.a, a)
  current.logf.b <- sigma.ab.fc.log(sigma.b, b)
  # Perform the jumps
  for (i in 1:iter) {
    # Get trial points
    new.sigma.a <- jump(sigma.a)
    new.sigma.b <- jump(sigma.b)
    # Get acceptance probabilities
    new.logf.a <- sigma.ab.fc.log(new.sigma.a, a)
    new.logf.b <- sigma.ab.fc.log(new.sigma.b, b)
    alpha.a <- exp(new.logf.a - current.logf.a)
    alpha.b <- exp(new.logf.b - current.logf.b)
    # Check acceptance
    if (runif(n=1, min=0, max=1) < alpha.a) {
      sigma.a <- new.sigma.a
      current.logf.a <- new.logf.a
    }
    if (runif(n=1, min=0, max=1) < alpha.b) {
      sigma.b <- new.sigma.b
      current.logf.b <- new.logf.b
    }
  }
  # Create a new matrix and return it
  return( matrix(c(sigma.a^2, 0, 0, sigma.b^2), ncol=2) )
}