aMTM: Adaptive Multiple-Try Metropolis algorithm
In fontaine618/aMTM: Adaptive Multiple-Try Metropolis Algorithm

Description Usage Arguments Details Value Author(s) References See Also Examples

View source: R/aMTM.R

This function performs the Adaptive Multiple-Try Metropolis algorithm as described in Fontaine and Bedard (2019). The sampling step is performed via a MTM algorithm with K gaussian proposals which may be correlated (see argument proposal and details) and using either weights that are proportional to the target density or importance weights. The adaptation step is performed via one of AM, ASWAM or RAM updates of the selected proposal density; a fixed global component may also be adapted at each iteration (see argument global and details) and the scale paramters may be adapted at each iteration (see argument scale and details). The AM and ASWAM update may be done using local steps rather than global steps (see argument local and details).

1	aMTM(target, N, K, x0, ...)

`target`	Target log-density which must be vectorized, i.e. take input of dimension Kd*. Additional parameters must be passed as a list trough the `parms` argument.
`N`	Size of the Monte Carlo sample.
`K`	Number of proposals in the MTM sampling.
`x0`	A vector of dimension d corresponding to the initial state of the chain.
`sig0`	An array of dimension ddK containing the K initial covariance of the instrumental gaussian distributions. Default is K identity matrices.
`mu0`	A matrix of dimension dK* containing the K initial mean parameters for AM and ASWAM updates. Default is K zero vectors.
`lam0`	A vector of dimension K containing the K scale parameters. Default is (2.38)^2/d for AM and AWSAM updates and 1 for RAM updates.
`adapt`	Determines the type of update of the covaraince done in the adaptation step : “AM” or 1 performs AM updates, “ASWAM” or 2 performs ASWAM updates, “RAM” or 3 performs RAM updates, any other value produces no update. Default is `"ASWAM"`.
`global`	Boolean value enabling the use of a global component that is not adapted (the first one). Default is `FALSE`.
`scale`	Boolean value enabling the adaptation of the scale parameter of the proposals that were not selected. It consists of reducing the scale aprameter whenever the running selection proportion is below 0.1/K. Default is `FALSE`.
`local`	Boolean value enabling the use of local update steps in AM or ASWAM updates. Default is `FALSE`.
`proposal`	Determines the type of proposals used in the MTM sampling : “ind” or 0 produces independant candidates, “common” or 1 produces candidates from a common random vectors, “QMC” or 2 produces candidates by a Randomized Quasi Monte Carlo procedure using a Koborov rule, “EA” or 3 produces extremely antithetic candidates. Default is `"ind"`.
`accrate`	Target acceptance rate for ASWAM and RAM updates. Default is 0.5.
`gamma`	Power used to produce the adaptation step size : at iteration `n`, the stepsize is given by `n^-gamma`. Default is 0.5 for RAM update and 0.7 otherwise and the range of suggested values is (0.5,1] to meet theoritical guarantees while values in (0,0.5] are accepted but may yield weird behaviour.
`parms`	A list of paramters passed to the `target` function for evaluation.
`weight`	Determines the type of weights used in the MTM sampling : “proportional” or 0 produces weights that are proportional to the target density, “importance” or -1 produces importance weights. Default is `proportional`.
`burnin`	A value in [0,1) indicating the amount of burnin to be done. The Monte Carlo sample size is adjusted to yield a sample of size N after the burn-in. Default is 0.

The MTM sampling with correlated candidates is described in Craiu and Lemieux (2007). The proposals all have Gaussian marginal distributions corresponding to a random walk (i.e. centered at the current state). The covariance is given by lambda^(k)Sigma^(k) where lambda^(k) is a positive scale parameter and Sigma^(k) is a symmetric positive-definite matrix.

The correlation between the proposals are defined by one of the four following procedures. First, the candidates may be independant (proposal="ind" or proposal=0) so that the joint density of the proposal is the product of the marginals. Second, the candidates may be construted from a common random vector (proposal="common" or proposal=1) in which case a d-dimensional standard noraml vector is generated and the candidates are all calculated from this random vector via the natural transformation. Third, the candidates may be generated via a randomized Qausi-Monte Carlo procedure using a Koborov rule (proposal="QMC" or proposal=2) meaning that they are regularly spaced, see Fontaine and Bedard (2019) for details. Fourth, the candidates may be extremely antithetic (proposal="EA" or proposal=3) which means that the correlation is chosen to maximize the Euclidian distance between the candidates, see Fontaine and Bedard (2019) for details.

The weight function used in the selection and acceptation steps may be one of the following two choices. First, the weight can be chosen to be proportional to the target density (weight="proportional" or weight=0), i.e. w^(k)(y|x)=pi(y). Second, the weights may take the form of importance weights (weight="imporatnce" or weight=-1), i.e.

w^(k)(y|x)=pi(y)/Q^(k)(y|x).

The adaptation of the parameters has different options. Unless no adaptation is performed, the update of the parameters of selected proposal density is performed at every iteration. If global=TRUE, then the first proposal density is updated at every iteration. If scale=TRUE, then all the scale parameters are adapted at every iterations to prevent that some of the proposals are never used. Now, the update of the parameters of a proposal may take one of three forms. First, the covariance Sigma^(k) can be adapted by the AM update of Haario et al. (2001) while the scale parameter remains constant (adapt="AM" or adapt=1). Second, the scale parameter may be adapted to attain a target acceptance rate as described by the ASWAM algorithm of Andrieu and Thoms (2008) where the covariance is adapted as in the AM update and the scale parameter is updated in the log scale to close in the gap between the mean accaptance rate and the target rate (adapt="ASWAM" or adapt=2). Third, the covariance may be updated by the RAM update of Vihola et al. (2012) where a target acceptance rate is used to adapt the covariance directly without updating the scale parameter (adapt="RAM" or adapt=3). In the AM and ASWAM update cases, it is possible to used local steps (local=TRUE) which may help the algorithm to adjust more precisely to local strutures of the target distribution. In particular, when the local moves are not used, all the covariances will often converge to the overall covariance of the target distribution.

In the context of adaptive MCMC using stochatic approximations to update the parameters, the adaptation stepsize gamma_n must satisfy some conditions. The algorithm uses gamma_n=n^-gamma for some parameter gamma. To insure good behaviour of the algorithm, we require 0<gamma<=1, but 0.5<gamma<=1 furthur garantees the convergence of the parameters so that the transition stabilizes in the long run.

The sampling and adaptation is done by a C++ core function to which this function is only a wrapper. The adaptation procedure rely on a rank one Cholesky update and downdate function taken from the RAMCMC package available at http://github.com/helske/ramcmc.

A list containing the following elements:

`X`	The MCMC sample in a matrix of dimension Nd* and is of class `"mcmc"` to insure compatibility with the `coda` package.
`acc.rate`	The acceptance rate of the sample.
`sel.prop`	The selection proportion of each proposal.
`mu`	The final value of the mean parameters (only meaningful for AM and ASWAM updates without local steps).
`lam`	The final value of the scale parameters.
`Sig`	The final value of the covariance paramters.
`sel`	A vector of length `N` containing the index of the selected proposal at each iteration.
`time`	The computing time.

Simon Fontaine, simfont@umich.edu

Andrieu, C. and Thoms, J. (2008). "A tutorial on adaptive MCMC". Statistics and computing, 18:4, 343-373.

Craiu, R.V. and Lemieux., C. (2007). "Acceleration of the multiple-try Metropolis algorithm using antithetic and stratifed sampling". Statistics and computing, 17:2, 109.

Fontaine, S. and Bedard, M. (2020+). "An Adaptive Multiple-Try Metropolis algorithm". To be submitted.

Haario, H., Saksman, E., Tamminen, J. et al. (2001). "An adaptive Metropolis algorithm". Bernoulli, 7:2, 223-242.

Helske, J. (2018). "ramcmc: Robust Adaptive Metropolis Algorithm". R package version 0.1.0-1, http://github.com/helske/ramcmc.

Liu, J.S., Liang, F. and Wong, W.H. (2000). "The Multiple-Try Method and Local Optimization in Metropolis Sampling". Journal of the American Statistical Association, 95:449, 121-134.

Vihola, M. (2012). "Robust adaptive Metropolis algorithm with coerced acceptance rate". Statistics and Computing, 22:5, 997-1008.

plot.aMTM for plotting and stats.aMTM for summaries.

## Not run: 
library(aMTM)
# Banana log-density with parameter B and a
p <- function(x, p) apply(x,1,function(x) -x[1]^2/(2*p$a^2) - 1/2*(x[2]+p$B*x[1]^2-p$B*p$a^2)^2)
# setup
set.seed(1)
N<-1e5;K<-3
B<-0.04;a<-8
# aMTM sampling with ASWAM update
mcmc <- aMTM(target=p, N=N, K=K, x0=c(0,0), parms=list(a=a,B=B), burnin=0.1)

# acceptance rate (target is 0.5)
mcmc$acc.rate

# plot 1D samples
plot(mcmc$X)
# plot 2D sample
plot(as.matrix(mcmc$X))
# add the final covariances
for(k in seq(K)) mixtools::ellipse(mcmc$mu[,k], mcmc$Sig[,,k]*mcmc$lam[k],alpha = 0.1, col=k)

# estimate of the mean (true value is (0,0))
colMeans(mcmc$X)
# estimate of the variance and true value
round(cov(mcmc$X),4)
matrix(c(a^2,0,0,1+B^2*a^4*2),2,2,T)

# estimation of the mean with MC standard error
mcmcse::mcse.mat(mcmc$X)
# asymptotic covariance of the estimator
mcmcse::mcse.initseq(mcmc$X)$cov
# multivariate effective sample size
mcmcse::multiESS(mcmc$X)

## End(Not run)