# R/metropdir.adapt.R In dirmcmc: Directional Metropolis Hastings Algorithm

#### Documented in metropdir.adapt

```#' @title Directional Metropolis Hastings with Adaptation.
#' @author Abhirup Mallik, \email{malli066@umn.edu}
#' @inheritParams metropdir
#' @param batchlen length of batch used for update. Default is 100.
#' @param targetacc Target acceptance ratio
#' @description Implements adaptive version of directional Metropolis Hastings.
#'@details This function is for automatically select a scaling factor for the directional Metropolis Hastings algorithm. This function uses batch wise update of the scale parameter to produce a adaptive chain. The user is required to supply a target acceptance ratio. The adaptive scheme modifies the scale parameter to achieve the target acceptance ratio. It is recommended that to check the output of the adaptation history as returned by the function.
#'
#' @examples
#' \dontrun{
#' Sigma <- matrix(c(1,0.2,0.2,1),2,2)
#' mu <- c(1,1)
#
#' Sig.Inv <- solve(Sigma)
#' Sig.det.sqrt <- sqrt(det(Sigma))
#
#' logf <- function(x,mu,Sig.Inv){
#'   x.center <- as.numeric((x-mu))
#'   out <- crossprod(x.center,Sig.Inv)
#'   out <- sum(out*x.center)
#'   -out/2
#'   }
#'
#' gr.logf <- function(x,mu,Sig.Inv){
#'   x.center <- as.numeric((x-mu))
#'   out <- crossprod(x.center,Sig.Inv)
#'   -as.numeric(out)
#' }
#' set.seed(1234)
#' system.time(out <- metropdir.adapt(obj = logf, dobj = gr.logf, initial = c(1,1),
#'                          lchain = 1e4,sd.prop=1,steplen=0,s=1, mu = mu,
#'                          Sig.Inv = Sig.Inv,targetacc=0.44))
#' #acceptance rate
#' out\$accept
#' #density plot
#' plot(density(out\$batch[,1]))
#' }
#' @seealso \code{\link{metropdir}} for DMH, \code{\link{iact}} for integrated auto correlation times, \code{\link{mcmcdiag}}, \code{\link{msjd}} for mean squared jump distance.
#' for summary of diagnostic measures of a chain, \code{\link{multiESS}} for Multivariate effective sample size.
#' @keywords dmh, mcmc, adaptive
#' @export

batchlen = 100,targetacc=0.234,...){
tol <- 1e-4
#sd.prop <- 1
p <- length(initial)
chain <- matrix(NA, nrow = lchain, ncol = p)
chain[1,] <- initial
accept <- 0

sdist2 <- function(x,s){
p <- length(x)
S <- diag(c(1/s,rep(1,p-1)))
crossprod(x,S)%*%x
}

x <- initial
dx <- dobj(x,...)
d.norm <- sqrt(sum(dx*dx))
if(d.norm > tol){
d <- dx/d.norm
}else d <- c(1,rep(0,p-1))
Gx <- qr.qy(qr(d),diag(p))
LAM <- diag(c(sqrt(s),rep(1,p-1)))
Sigma.sqrt.x <- Gx%*%LAM%*%t(Gx)
mu.x <- x+d*rep(steplen,p)*d.norm
acceptance <- 0

#Initialize for adaption #
bc <- 1
bi <- 1
bataccept <- 0
numbatch <- floor(lchain/batchlen)
sdout <- numeric(numbatch)

for(i in 1:lchain){
chain[i,] <- x
z <- rnorm(p,sd=sd.prop)
mu.x <- x + dx*rep(steplen,p)
y <- Sigma.sqrt.x%*%z+mu.x

dy <- dobj(y,...)
d.norm <- sqrt(sum(dy*dy))
if(d.norm > tol){
d <- dy/d.norm
}else d <- c(1,rep(0,p-1))
Gy <- qr.qy(qr(d),diag(p))
mu.y <- y+d*rep(steplen,p)*d.norm
Sigma.sqrt.y <- Gy%*%LAM%*%t(Gy)

##diff <- x-y
logalpha <- obj(y,...) - obj(x,...) +
(sdist2(Gx%*%(y-mu.x),s)-sdist2(Gy%*%(x-mu.y),s))/2

if(logalpha > 0){
accept <- TRUE
} else {
u <- runif(1)
if(logalpha > log(u)){
accept <- TRUE
} else {
accept <- FALSE
}
}

if(accept){
x <- y
dx <- dy
mu.x <- mu.y
Gx <- Gy
Sigma.sqrt.x <- Sigma.sqrt.y
acceptance <- acceptance + 1
}
if(bi < batchlen){
if(accept) bataccept <- bataccept + 1
}else{
sdout[bc] <- sd.prop
if(accept) bataccept <- bataccept + 1
if(bataccept > batchlen*targetacc){
sd.prop <- sd.prop * exp(min(c(0.01,1/sqrt(bc))))
} else {
sd.prop <- sd.prop / exp(min(c(0.01,1/sqrt(bc))))
}
bi <- 0
bc <- bc + 1
bataccept <- 0
}
bi <- bi + 1
}
list(accept = acceptance/lchain, batch = chain,sdout = sdout)
}
```

## Try the dirmcmc package in your browser

Any scripts or data that you put into this service are public.

dirmcmc documentation built on May 2, 2019, 4:14 a.m.