metropdir {dirmcmc}

Directional Metropolis Hastings

Description

Implements Markov Chain Monte Carlo for continuous random vectors using Directional Metropolis Hasting Algorithm

Usage

metropdir(obj, dobj, initial, lchain, sd.prop = 1, steplen = 0, s = 0.95,
  ...)

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

...

any arguments to be passed to obj and dobj.

Details

Runs a “Directional Metropolis Hastings” algorithm, with multivariate normal proposal producing a Markov chain with equilibrium distribution having a specified unnormalized density. Distribution must be continuous. Support of the distribution is the support of the density specified by argument obj.

Value

Returns the following objects in a list:

• accept. acceptance rate.

• batch. resulting chain.

Author(s)

Abhirup Mallik, [e-mail]

See Also

metropdir.adapt for adapting 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(obj = logf, dobj = gr.logf, initial = c(1,1),
                         lchain = 1e4,sd.prop=1,steplen=0,s=1, mu = mu,
                         Sig.Inv = Sig.Inv))
#acceptance rate
out$accept
#density plot
plot(density(out$batch[,1]))

## End(Not run)