atmcmc: Runs a MCMC algorithm tuned via adaption and corresponding...

Description Usage Arguments Details Value References Examples

Description

A symmetric random walk Metropolis algorithm with Gaussian proposals is automatically tuned and run and diagnosed to converge to a target distribution and estimate the stationary mean of a functional

Usage

1
2
3
4
5
6
7
8
atmcmc(g, dim, X0, support = cbind(rep(-Inf, dim), rep(Inf, dim)),
  multimodal = F, functional = function(X) {     X }, maxiter = 2e+06,
  mult = (2.38)^2/(dim), mrep = 10, nrep = 10, batchwidth = 200,
  holdup = 10, batchwidth.adp1 = 100, scale.adj = 0.05,
  endbatch.adp1 = 2, minaccpt = 0.28, maxaccpt = 0.6, nreg = 5,
  startdist = 1.5, minR = 0.9, maxR = 1.1, CI.alpha = 0.05,
  nimprob.X = 100, minaccpt.adp2 = 0.02, batchwidth.adp2 = 200,
  jumpprob = 0.05, displayfreq = 100, plot = T, m = 10)

Arguments

g

log of target density function

dim

dimension of target density function

X0

initial value of Markov chain

support

support of target density function. Argument takes a ‘dim’ x 2 matrix

multimodal

whether to assume target density function is strongly multimodal. Argument takes T or F

functional

function to estimate an expectation with respect to target density function

maxiter

maximum iteration number for a full MCMC run

mult

constant multiplied to the covariance matrix of the proposal distribution in the 2nd adaption phase and the sampling phase

mrep

number of chains created to find multiple modes if multimodal=T

nrep

number of replicative chains for Gelman-Rubin diagnostics (convergence diagnostics)

batchwidth

number of iterations in a batch. Is used to find the average X values in the transient phase and the average squared jumping distances in the 2nd adaption phase.

holdup

holdup*batchwidth is the number of iterations between the start of the sampling phase and the first computation of Gelman-Rubin diagnostics (to check for convergence of the chain) in the sampling phase

batchwidth.adp1

number of iterations in a batch initially used to check for the acceptance rate in the 1st adaption phase

scale.adj

amount added to or subtracted from the log of standard deviation of the univariate Gaussian proposal for each coordinate in the 1st adaption phase

endbatch.adp1

batchwidth.adp1 x 2^(endbatch.adp1) iterations needs to have the desired level of acceptance rate to end the 1st adaption phase

minaccpt

minimum acceptable acceptance rate to end the 1st adaption phase

maxaccpt

maximum acceptable acceptance rate to end the 1st adaption phase

nreg

number of distinct points in the simple linear regression to check if there is a linear trend in average X values in the transient phase or in average squared jumping distances in the 2nd adaption phase

startdist

number multiplied to ‘max X value - min X value’ in the 2nd adaption + flat part of transient phase (or just 2nd adaption phase if multimodal=T). Is used to determine over-disposed starting distribution for the sampling phase

minR

minimum acceptable R value to end the sampling phase

maxR

maximum acceptable R value to end the sampling phase

CI.alpha

(1-CI.alpha)x100% confidence interval is constructed for R_interval in the sampling phase

nimprob.X

number of consecutive g(X)=-Inf iterations required to break off the full algorithm. This is to prevent the chain from drifting to a wrong direction

minaccpt.adp2

minimum acceptable acceptance rate in the 2nd adaption phase to keep the value of ‘mult’ as it is. If the acceptance rate is less than ‘minaccpt.adp2’, ‘mult’ is cut down to 'mult'/max(2,'dim')

batchwidth.adp2

first n number of iterations used in the 2nd adaption phase to check for the acceptance rate to decide whether to decrease ‘mult’

jumpprob

probability of a jump between modes at each iteration in the sampling phase

displayfreq

how frequently to display the iteration number as ‘atmcmc’ runs. Every ‘displayfreq’th iteration number is printed on the screen

plot

whether to display traceplots of coordinate 1 (& mode 1 for multimodal=T) as each phase ends. Takes argument T or F

m

every ‘m’th iteration is plotted. Has to be a factor of ‘batchwidth’/2

Details

The algorithm automatically tunes the covariance matrix of the proposal distribution N(X_n, Σ_p) of a symmetric random walk Metropolis algorithm. The algorithm can be broken down into four main phases: a 1st adaption phase, transient phase, 2nd adaption phase and sampling phase. The 1st adaption phase employs the Adaptive Metropolis-within-Gibbs algorithm from Roberts and Rosenthal (2009), and the diagnostics to end this phase is to check whether the acceptance rate for every coordinate comes into the desired range. The transient phase runs a Metropolis-within-Gibbs algorithm, and this runs until the chain reaches the mode of the target distribution. The purpose of the transient phase is to avoid including bad X values when tuning for Σ_p = mult * Σ_n, where Σ_n is the empirical covariance matrix of the target distribution calculated from the values generated by the Markov chain. The diagnostics to end this phase is to fit a simple linear regression to see if the chain values are trending. The 2nd adaption phase employs a slightly modified version of the Adaptive Metropolis algorithm from Haario et al. (2001) or Roberts and Rosenthal (2009). This phase updates Σ_p at every iteration by calculating the empirical covariance matrix of the target distribution from the chain values. The diagnostics to confirm whether this phase is indeed improving the chain is to see if the squared jumping distance between every consecutive iteration is increasing. Again, a simple linear regression is used to see if the squared jumping distances are increasing. After all this, a symmetric random walk Metropolis algorithm is run and Gelman-Rubin diagnostics is used to verify convergence of the Markov chain. Note that 2nd half of the sampling phase is what we take as a sample.

For a target distribution that is considered to be ‘strongly multimodal’, the basic structure of the algorithm is still the same, but multiple chains are run in the 1st adaption phase and transient phase until each chain reaches different mode. The algorithm leaves only one chain for each unique local mode and deletes others. It considers two modes are different when, in at least one coordinate, the absolute value of the difference of two means is lesser than the smaller of the standard deviation of two. A 2nd adaption phase is run for each remaining chain, and after the 2nd adaption phase, the algorithm confirms whether each chain has different mode. For the sampling phase, at each iteration, the chain either moves inside one local mode or jumps to another mode at a fixed probability. Again, Gelman-Rubin diagnostics is used to check for convergence.

Value

A list consisting of

Xveclistdim1, Xveclistdim2,...

matrix of X values saved from the sampling phase. Each matrix contains X values of each coordinate

Xveclistbase

matrix of X values saved from the 1st adaption ,transient, and 2nd adaption phase. Includes only one chain for each unique local mode for multimodal=T

nummodes

number of unique local modes found for multimodal=T

dim

dimension of target density function

batchwidth

number of iterations in a batch. Is used to find the average X values in the transient phase and the average squared jumping distances in the 2nd adaption phase.

means

average of X values from the 2nd half of sampling phase

functionalmeans

average of functional X values from the 2nd half of sampling phase

sumchain

string of iteration numbers to show when each phase has ended. For the sampling phase, this shows when the 1st half of sampling phase has ended also

acceptrate

acceptance rate of the 2nd half of sampling phase

runtime

runtime of the full MCMC

It also prints values of estimates, estimates_of_functional, acceptance_rate, time_elapsed, and phase_length. For details, see ‘summarymcmc’ section

References

Heikki Haario, Eero Saksman, and Johanna Tamminen. An adaptive metropolis algorithm. Bernoulli, 7(2):223-242, 2001.

Gareth O. Roberts and Jeffrey S. Rosenthal. Examples of adaptive MCMC. Journal of Computational and Graphical Statistics, 18(2):349-367, 2009.

Andrew Gelman and Donald B. Rubin. Inference from iterative simulation using multiple sequences. Statistical science, 7(4):457-472, 1992.

Stephen P. Brooks and Andrew Gelman. General methods for monitoring convergence of iterative simulations. Journal of computational and graphical statistics, 7(4):434-455, 1998.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
dim = 3  #dimension of target density function
X0 = rep(0.1,dim)  #initial X value

tmpmat = rbind(c(-0.7, 1.2, 1.6),c(0.9, 1.1, -0.1),c(0.2, 0.3, -1.5))
targSigma = t(tmpmat) %*% tmpmat
targMean = c(22, -10,  15)
#log of target density function
g = function(X){-0.5*t(X-targMean)%*%solve(targSigma)%*%(X-targMean)}

output = atmcmc(g, dim, X0)
plotmcmc(output, name = "project1")
summarymcmc(output, name = "project1")

atmcmc documentation built on May 2, 2019, 11:05 a.m.