Nothing
### A very quick demonstration of no-U-turn sampling in ADMB
library(adnuts)
### ----------- ADMB example
## This is the packaged simple regression model
path.simple <- system.file('examples', 'simple', package='adnuts')
## It is best to have your ADMB files in a separate folder and provide that
## path, so make a copy of the model folder locally.
path <- 'simple'
dir.create(path)
trash <- file.copy(from=list.files(path.simple, full.names=TRUE), to=path)
## Compile and run model
setwd(path)
## If admb is not in the PATH this line will fail and you need to manually compile the model.
system('admb simple.tpl')
system('simple')
setwd('..')
## 3 options to specify inits: list, NULL (uses MLE), or a function that
## returns a list.
init <- lapply(1:3, function(i) rnorm(2))
init <- NULL # uses MLEs -- not recommended!
init <- function() rnorm(2)
fit <- sample_nuts(model='simple', init=init, path=path, cores=1)
fit$cmd[1] # this is the command line arguments used
## Merged chains after discarding warmup phase
post <- extract_samples(fit)
str(post)
## A list with MLE fit
str(fit$mle)
## Can also run in parallel which is the default when
## chains>1. Here we also execute the mceval phase
## afterward. Note that the .psv file only has post-warmup and
## thinned (if used) samples in it. Chains are rbind'ed together,
## and thus the -mceval call runs on all chains and can be used
## directly for inference. You can set mceval=TRUE here, or run
## your model later manually as typically done.
fit <- sample_nuts(model='simple', init=init, path=path, mceval=TRUE)
## Can also specify a duration argument for capping the run
## time. This is useful e.g. if you want to run overnight but
## have results by 8am. Here we do 0.5 minutes just to
## demonstrate. This period needs to be long enough to do the
## warmup or it'll throw an error. Normally we wouldn't use
## thinning for NUTS but only to demonstrate here since the model
## runs so fast. I recommended setting a reasonable warmup, then
## an unreasonably big iter. It will truncate to the shortest
## chains.
fit <- sample_nuts(model='simple', init=init, path=path,
warmup=200, iter=20000000, thin=100,
duration=.5)
## The default is use an adaptive mass matrix. If we want to use
## the MLE covariance as the mass matrix, set it using
## control. This *should* make for more efficient sampling in
## most cases. For technical reasons, must optimize the model
## with flag -hbf 1. Here setting skip_optimization=FALSE will
## rerun the model before starting the chains. Or you could
## manually optimize the model with the flag '-hbf'.
fit <- sample_nuts(model='simple', init=init, path=path,
skip_optimization=FALSE,
control=list(metric='mle'))
## See the vignette for a discussion on the options for the
## metric (mass matrix adaptation and such).
## Can also use a slightly modified version of the original
## Metropolis MCMC algorithm (Random Walk Metropolis or
## RWM). Here we really want to use the MLE covariance. Since we
## ran above with hbf=1 we need to rerun model to recreate the
## .hes and .cov files. Can do this with skip_optimization
## argument again.
fit.rwm <- sample_rwm(model='simple', init=init, path=path,
skip_optimization=FALSE,
iter=200000, thin=100, mceval=TRUE)
### ------------------------------------------------------------
### Convergence diagnostics
## Key convergence diagnostics (effective sample size and Rhat)
## are available in the summary of the fit. It is your
## responsibility to check for signs of non-convergence before
## using the output for inference!!
summary(fit)
## or from the rstan::monitor output
str(fit$monitor)
### ------------------------------------------------------------
### Plotting options
pairs_admb(fit) # modified pairs just for ADMB fits like this
## Can also use ShinyStan (make sure to exit it)
## launch_shinyadmb(fit)
plot_sampler_params(fit) # NUTS adaptation
## Compare MLE and posterior of marginals. See help for
## recommendation for creating multipage PDF for high dimensional
## parameters.
plot_marginals(fit)
### ------------------------------------------------------------
### Extracting posterior samples
## Get post-warmup, merged chains into data frame (these two are
## identical)
str(as.data.frame(fit))
str(extract_samples(fit))
## If you want it in list form, e.g., to put into coda package
str(extract_samples(fit, as.list=TRUE))
## If you want to see warmup samples, and the log-posterior (lp__ column)
str(extract_samples(fit, inc_warmup=TRUE, inc_lp=TRUE))
## Remove folder
unlink(path, TRUE)
### End of demo
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.