adapt_S: Update the Proposal of RAM Algorithm

In helske/ramcmc: Robust Adaptive Metropolis Algorithm

Description

Given the lower triangular matrix S obtained from the Cholesky decomposition of the shape of the proposal distribution, function adapt_S updates S according to the RAM algorithm.

Usage

adapt_S(S, u, current, n, target = 0.234, gamma = 2/3)

Arguments

S	A lower triangular matrix corresponding to the Cholesky decomposition of the scale of the proposal distribution.
u	A vector with with length matching with the dimensions of S.
current	The current acceptance probability.
n	Scaling parameter corresponding to the current iteration number.
target	The target acceptance rate. Default is 0.234.
gamma	Scaling parameter. Default is 2/3.

Value

If the resulting matrix is positive definite, an updated value of S. Otherwise original S is returned.

Note

If the downdating would result non-positive definite matrix, no adaptation is performed.

References

Matti Vihola (2012). "Robust adaptive Metropolis algorithm with coerced acceptance rate". Statistics and Computing, 22: 997. doi:10.1007/s11222-011-9269-5

Examples

# sample from standard normal distribution
# use proposals from the uniform distribution on
# interval (-s, s), where we adapt s

adapt_mcmc <- function(n = 10000, s) {
  x <- numeric(n)
  loglik_old <- dnorm(x[1], log = TRUE)
  for (i in 2:n) {
    u <- s * runif(1, -1, 1)
    prop <- x[i] + u
    loglik <- dnorm(prop, log = TRUE)
    accept_prob <- min(1, exp(loglik - loglik_old))
    if (runif(1) < accept_prob) {
      x[i] <- prop
      loglik_old <- loglik
    } else {
      x[i] <- x[i - 1]
    }
    # Adapt only during the burn-in
    if (i < n/2) {
      s <- adapt_S(s, u, accept_prob, i)
    }
  }
  list(x = x[(n/2):n], s = s)
}

out <- adapt_mcmc(1e5, 2)
out$s
hist(out$x)
# acceptance rate:
1 / mean(rle(out$x)$lengths)

helske/ramcmc documentation built on Oct. 9, 2021, 9:38 a.m.