Implementation of the Adaptive Gibbs and Adaptive Metropolis within Adaptive Gibbs Samplers described in C., Latuszynski, and Roberts (2018).
The package provides a set of functions to run the Adaptive Random Scan Gibbs (ARSG) and Adaptive Metropolis-within-Gibbs (AMwG) algorithms for a user-defined target distribution. The adaptive algorithms are defined in the ARSG paper. The adaptations tune the sampling weights of the algorithms, and also the variances of the proposals in the AMwG algorithm. The algorithms are implemented in C++ providing R interface. While the user can define the target density in R, the performance of the algorithm is much better if the density is described in C++ using the provided template. Gibbs Sampling is implemented only for C++ defined densities, where the user has to explicitly define a function that samples from full conditionals. The main function AMCMC(...) produces a trace of the chain which is saved as a set of text files in a pre-specified folder. Please refer to appendix.txt for details for additional functions provided by the package.
This package was tested on OS X 10.10 and Ubuntu 16.04
Mac Users:
- Make sure the newest version of R is installed. If RStudio is used, please update to at least version 0.10.3
- Open terminal:
- Run
xcode-select --install
to install xcode command line tools (this will install gcc in particular). - If gfortran is not installed on your machine, run
curl -O http://r.research.att.com/libs/gfortran-4.8.2-darwin13.tar.bz2
sudo tar fvxz gfortran-4.8.2-darwin13.tar.bz2 -C /
This should update your gfortran compiler (see here for details)
Ubuntu Users:
- Make sure the last version of R is installed.
- Make sure you have installed gcc:
sudo apt-get install gcc
Install Rcpp, RcppArmadillo, RcppParallel, if not yet installed:
- Open R/RStudio session and install the following packages if not yet installed:
install.packages("Rcpp","RcppArmadillo", "RcppParallel")
-
Download/clone the repository files to a local directory, say, "~/Downloads".
-
Open R/Rstudio session use
Rcpp::sourceCpp
to compile the library:
Rcpp::sourceCpp("~/Downloads/Adaptive-Gibbs-Sampler/Adaptive_Gibbs.cpp")
- If you want to compile the library later to a specific folder, say, "~/Downloads/Adaptive-Gibbs-SamplerAdaptiveGibbsLib", you need to specify
cacheDir
argument:
Rcpp::sourceCpp("~/Downloads/Adaptive-Gibbs-Sampler/Adaptive_Gibbs.cpp", cacheDir = "~/Downloads/Adaptive-Gibbs-Sampler/AdaptiveGibbsLib")
- Use Step 3 to load precompiled library. It will be recompiled automatically if any changes are made to the source files.
After compiling the library:
- First, set up a folder where all the output will be stored using
set_working_directory(path)
, wherepath
is the desired folder, say,
set_working_directory("Users/Admin/AdaptiveGibbs/simulation_results/")
Notice "/" at the end of the pass.
- Now you can use
AMCMC(...)
and related functions. Please refer to AMCMC_info.md, AMCMC_appendix.md, and tutorials to learn how to useAMCMC(...)
function.
Consider sampling from a 10-dimensional multivariate normal target distribution with a covariance matrix S
. As an example, the user can generate a block diagonal matrix using set_example_covariance
:
dim <- 10 # dimensionality of the target
N <- 10000 # number of desired samples
# set an example covariance matrix
set_example_covariance(dim)
# print the matrix on screen
get_covariance()
# set random seed
set.seed(42)
In order to sample from the distribution described in C++ language, at least its reference name should be provided. The gaussian target distribution is described in gaussian_target.hpp file. Sampling can be performed using
adaptive_chain<-AMCMC(distribution_type = "gaussian", logdensity = 1, dimension = dim, N = N)
The function calls the coordinate-wise Adaptive Metropolis within Adaptive Gibbs Sampler for gaussian distribution. logdensity
indicates that the log-density of the target distribution is used for Metropolis acceptance-rejection step. The output of the chain is written to the prespecified directory (use get_working_directory()
to print the directory on screen). adaptive_chain
stores the value of estimated inverse pseudo-spectral gap and the optimal sampling selection probabilities.
Use trace_coord(i)
to access the trace of the coordinate i
. For example, in order to estimate the covariance matrix, one could use the following code
S <- matrix(0,nrow = N, ncol = dim)
for(i in 1:dim)
{
v <- trace_coord(i)
S[,i] <- v
S <- cor(S)
# estimated correlation matrix:
S
Detailed guide on writing user-defined densities is described in write_custom_density.md.
In short, one could provide an R-defined density, say,
example_logdensity<-function(x)
{
return( -1./2* (t(x) %*%Q_0%*% x)[1,1] )
}
and call
AMCMC(R_density = example_logdensity, logdensity = 1, dimension = dim, N = N)
However, a much faster performance is achieved if a C++ density is supplied. A template file is provided in template.hpp. You have to specify at least one function of the template class, say, log-density function:
double gaussian::logdensity(vec theta)
{
return -1./2*dot(theta.t()*Q,theta);
}
Precision matrix Q
should be specified in the constructor.
For various examples follow this link.