metropdir: Directional Metropolis Hastings

In dirmcmc: Directional Metropolis Hastings Algorithm

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)

Example output

   user  system elapsed 
  1.170   0.004   1.187 
[1] 0.5405

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