metropdir.adapt {dirmcmc}

Directional Metropolis Hastings with Adaptation.

Description

Implements adaptive version of directional Metropolis Hastings.

Usage

metropdir.adapt(obj, dobj, initial, lchain, sd.prop = 1, steplen = 0,
  s = 0.95, batchlen = 100, targetacc = 0.234, ...)

Arguments

obj

an R function that evaluates the log unnormalized probability density of the desired equilibrium distribution of the Markov chain. First argument is the state vector of the Markov chain. Other arguments arbitrary and taken from the ... arguments of this function. Should return - Inf for points of the state space having probability zero under the desired equilibrium distribution.

dobj

an R function that evaluates the derivative of the log unnormalized probability density at the current state of the markov chain.

initial

Initial state of the markov chain. obj(state) must not return #' -Inf

lchain

length of the chain

sd.prop

scale to use for the proposal

steplen

tuning parameter in mean of proposal

s

tuning parameter in the covariance of proposal

batchlen

length of batch used for update. Default is 100.

targetacc

Target acceptance ratio

...

any arguments to be passed to obj and dobj.

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.

Author(s)

Abhirup Mallik, [e-mail]

See Also

metropdir for DMH, iact for integrated auto correlation times, mcmcdiag, msjd for mean squared jump distance. for summary of diagnostic measures of a chain, multiESS for Multivariate effective sample size.

Examples

## Not run: 
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]))

## End(Not run)