knitr::opts_chunk$set(fig.path = "man/figures/", warning = FALSE)
The fmcmc
R package provides a lightweight general framework for implementing
Markov Chain Monte Carlo methods based on the Metropolis-Hastings algorithm. This
implementation's primary purpose lies in the fact that the user can incorporate the
following flexibly:
Automatic convergence checker: The algorithm splits the MCMC runs according
to the frequency with which it needs to check convergence. Users can
use either one of the included functions (convergence_gelman,
convergence_geweke,
etc.), or provide their own.
Run multiple chains in parallel fashion: Using either a PSOCK
cluster
(default), or providing a personalized cluster object like the ones in
the parallel
R package.
User defined transition kernels: Besides of canonical Gaussian Kernel,
users can specify their own or use one of the included in the package, for
example: kernel_adapt
, kernel_ram
, kernel_normal_reflective
, kernel_unif
,
kernel_mirror
, or kernel_unif_reflective
.
All the above without requiring compiled code. For the latest about fmcmc,
checkout
the NEWS.md section.
While a lot of users rely on MCMC tools such as Stan (via the rstan package) or WinBUGS (via rstan), in several settings either these tools are not enough or provide too much for things that do not need that much. So, this tool is for you if:
You have a simple model to estimate with Metropolis-Hastings.
You want to run multiple chains of your model using out-of-the-box parallel computing.
You don't want (or cannot) rely on external tools (so you need good-old base R only for your models).
You want to implement a model in which your model parameters are either bounded (like a standard error, for example), or are not, say, continuous (e.g., a size variable in a Binomial distribution), so you need a personalized transition kernel.
In any other case, you may want to take a look at the previously mentioned R packages, or check out the mcmc R package, which also implements the Metropolis-Hastings algorithm (although with not all the features that this R package has), the adaptMCMC R package, or the MCMCpack R package.
If you want to get the latest bleeding-edge version from Github, you can use devtools:
devtools::install_github("USCbiostats/fmcmc")
The latest (stable) release is also available on CRAN:
install.packages("fmcmc")
citation("fmcmc")
In the following we show how to use the package for estimating parameters in a linear regression model. First, let's simulate some data to use:
set.seed(78845) n <- 1000 X <- rnorm(n) y <- 3.0 + 2.0*X + rnorm(n, sd = 4)
In this case, we have three parameters to estimate, the constant (2.0), the $\beta$ coefficient (2.0), and the standard deviation parameter of the error (1.5).
To estimate this model, we can either maximize the log-likelihood function--which is what is usually done--or we could do it using MCMC. In either case, we need to specify the log(unnormalized some times)-likelihood function:
ll <- function(p, X., y.) { joint_ll <- dnorm(y. - (p[1] + X.*p[2]), sd = p[3], log = TRUE) joint_ll <- sum(joint_ll) # If is undefined, then we explicitly return infinte (instead of NaN, for # example) if (!is.finite(joint_ll)) return(-Inf) joint_ll }
Notice that the function has more than one argument, in this case, p,
the vector of parameters, X.
and y.,
the data of our model.
Let's do a first run of the MCMC algorithm using the function of the same name (first, load the package, of course):
library(fmcmc) # Running the MCMC (we set the seed first) set.seed(1215) ans <- MCMC( ll, initial = c(0, 0, sd(y)), nsteps = 5000, X. = X, y. = y )
As the output object is an object of class mcmc
from the coda
R package, we
can use all the functions from it on our output:
library(coda) plot(ans) summary(ans)
While the summary statistics look very good (we got very close to the original
parameters), the trace of the parameters looks very bad (poor mixing). We can
re-run the algorithm changing the scale parameter in the kernel_normal
function.
To do so, we can pass ans
as the initial
argument so that
the function starts from the last point of that chain:
ans <- MCMC( ll, initial = ans, nsteps = 5000, X. = X, y. = y, kernel = kernel_normal(scale = .05) # We can set the scale parameter like this ) plot(ans)
Much better! Now, what if we use Vihola (2012) Robust Adaptive Metropolis (which is also implemented in the R package adaptMCMC)
ans_RAM <- MCMC( ll, initial = ans, nsteps = 5000, X. = X, y. = y, kernel = kernel_ram() ) plot(ans_RAM) 1 - rejectionRate(ans_RAM)
We can also try using Haario et al. (2001) Adaptive Metropolis
ans_AM <- MCMC( ll, initial = ans, nsteps = 5000, X. = X, y. = y, kernel = kernel_adapt() ) plot(ans_AM) 1 - rejectionRate(ans_AM)
Finally, if needed, we can also access information about the last run using MCMC_OUTPUT
.
For example, if we wanted to look at the trace of the logposterior function, we could
use the get_logpost()
function:
plot(get_logpost(), type = "l")
The set of proposed values is also available using the get_draws()
function:
boxplot(get_draws(), type = "l")
If the previous run featured multiple chains, then get_logpost()
would return a list
instead of length get_nchains()
.
Now, suppose that the algorithm actually takes a lot of time to actually reach stationary state, then it would be nice to actually sample from the posterior distribution once convergence has been reached. In the following example we use multiple chains and the Gelman-Rubin diagnostic to check for convergence of the chain:
set.seed(1215) # Same seed as before ans <- MCMC( ll, initial = c(0, 0, sd(y)), nsteps = 5000, X. = X, y. = y, kernel = kernel_normal(scale = .05), nchains = 2, # Multiple chains conv_checker = convergence_gelman(200) # Checking for conv. every 200 steps )
As a difference from the previous case, now we didn't have to wait until the 5,000 completed, but the algorithm stopped for us, allowing us to start generating the desired sample much quicker.
For this final example, we will use the kernel
argument and provide what
corresponds to a transition kernel which makes proposals within certain boundaries,
in particular, we want the algorithm to propose only positive values for the sd
parameter, which we now must be positive.
Moreover, since we know that we will only get positive values, we
can go further and modify ll
skipping the check for finite values:
ll <- function(p, X., y.) { sum(dnorm(y. - (p[1] + X.*p[2]), sd = p[3], log = TRUE)) }
Much simpler function! Let's do the call of the MCMC function specifying the
right transition kernel to increase the acceptance rate. In this example, we will
set the max of all parameters to be 5.0, and the min to be -5.0 for the constant
and 0 for the beta coefficient and the variance parameter, all this using the
kernel_normal_reflective
(which implements a normal kernel with boundaries) function:
set.seed(1215) # Same seed as before ans <- MCMC( ll, initial = c(0, 0, sd(y)), nsteps = 5000, X. = X, y. = y, kernel = kernel_normal_reflective( ub = 5.0, # All parameters have the same upper bound lb = c(-5.0, 0.0, 0.0), # But lower bound is specified per parameter scale = 0.05 # This is the same scale as before ), nchains = 2, conv_checker = convergence_gelman(200) )
Again, as the proposal kernel has lower and upper bounds, then we are guaranteed that all proposed states are within the support of the parameter space.
fmcmc
is just one way to work with Markov Chain Monte Carlo. Besides stan
and
WinBUGS
, there are other ways to do MCMC in R: mcmc,
HybridMC, adaptMCMC, and
elhmc among others (take a look at
the CRAN Task View on Bayesian Inference.)
fmcmc
We welcome contributions to fmcmc
. Whether reporting a bug, starting a discussion by asking a question, or proposing/requesting a new feature, please go by creating a new issue here so that we can talk about it.
Please note that the 'fmcmc' project is released with a Contributor Code of Conduct. By contributing to this project, you agree to abide by its terms.
Supported by National Cancer Institute Grant #1P01CA196596.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.