| Title: | Robust Adaptive Metropolis Algorithm |
|---|---|
| Description: | Function for adapting the shape of the random walk Metropolis proposal as specified by robust adaptive Metropolis algorithm by Vihola (2012) <doi:10.1007/s11222-011-9269-5>. The package also includes fast functions for rank-one Cholesky update and downdate. These functions can be used directly from R or the corresponding C++ header files can be easily linked to other R packages. |
| Authors: | Jouni Helske [aut, cre] |
| Maintainer: | Jouni Helske <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.1.2 |
| Built: | 2024-11-02 03:57:11 UTC |
| Source: | https://github.com/helske/ramcmc |
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.
adapt_S(S, u, current, n, target = 0.234, gamma = 2/3)adapt_S(S, u, current, n, target = 0.234, gamma = 2/3)
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. |
If the resulting matrix is positive definite, an updated value of S. Otherwise original S is returned.
If the downdating would result non-positive definite matrix, no adaptation is performed.
Matti Vihola (2012). "Robust adaptive Metropolis algorithm with coerced acceptance rate". Statistics and Computing, 22: 997. doi:10.1007/s11222-011-9269-5
# 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)# 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)
Given the lower triangular matrix L obtained from the Cholesky decomposition of A,
function chol_downdate updates L such that it corresponds to the decomposition of A - u*u'
(if such decomposition exists).
chol_downdate(L, u)chol_downdate(L, u)
L |
A lower triangular matrix. Strictly upper diagonal part is not referenced. |
u |
A vector with with length matching with the dimensions of L. |
Updated L.
The function does not check that the resulting matrix is positive semidefinite.
Given the lower triangular matrix L obtained from the Cholesky decomposition of A,
function chol_update updates L such that it corresponds to the decomposition of A + u*u'.
chol_update(L, u)chol_update(L, u)
L |
A lower triangular matrix. Strictly upper diagonal part is not referenced. |
u |
A vector with with length matching with the dimensions of L. |
Updated L.
L <- matrix(c(4,3,0,5), 2, 2) u <- c(1, 2) chol_update(L, u) t(chol(L %*% t(L) + u %*% t(u)))L <- matrix(c(4,3,0,5), 2, 2) u <- c(1, 2) chol_update(L, u) t(chol(L %*% t(L) + u %*% t(u)))