#'
#' Fit \code{bdm} model
#'
#' Execute a Bayesian model fit using \pkg{rstan}.
#'
#' By default a Bayesian fit is executed through a call to \code{\link[rstan:sampling]{sampling}}, which implements an MCMC algorithm. Default values for \code{chains}, \code{iter}, \code{warmup} and \code{thin} follow those for \pkg{rstan}.
#'
#' The \code{init} argument can be a \code{list}, \code{function} or \code{character} string. If it is a function then it should take no arguments and return a named list of intial values for the estimated parameters. Alternatively the list can be specified directly.
#' This behaviour matches that for \code{\link[rstan:sampling]{sampling}}. If a character string is supplied it should be either \code{'random'} or \code{'fixed'}.
#' If the model is the default model and \code{init = 'fixed'} then sensible starting values for \code{r}, \code{logK} and \code{x} are produced using \code{\link{getr}}, \code{\link{getlogK}} and \code{\link{getx}}.
#' If the model is the default model and \code{init = 'random'} then sensible starting values are obtained by sampling from the priors for \code{r, logK, x}.
#' If the model is not the default model, then the user should specify a function or list, otherwise starting values will be randomly generated by \code{\link[rstan:sampling]{sampling(init = 'random', ...)}}.
#'
#' @param object a \code{bdm} model object
#' @param data a \code{list} object containing the model inputs
#' @param run optional character vector to label the run
#' @param init an initialisation \code{list}, \code{function} or \code{character} string
#' @param chains number of MCMC chains
#' @param iter number of iterations per chain
#' @param warmup number of iterations to be discarded
#' @param thin sampling interval from chains
#' @param ... further arguments to \code{\link[rstan:sampling]{sampling}}
#'
#' @return Returns a \code{bdm} object containing posterior samples contained in \code{object@@trace}.
#'
#' @examples
#' # get some data
#' data(albio)
#' dat <- bdmData(harvest = albio$catch, index = albio$cpue, time = rownames(albio))
#'
#' # initialize and fit default model
#' \dontrun{
#' mdl <- bdm()
#' mdl <- compiler(mdl)
#' mdl <- sampler(mdl, dat)
#' }
#'
#' @include bdm.R
#'
#' @import methods
#' @import rstan
#'
#' @export
setGeneric("sampler", function(object, ...) standardGeneric("sampler"))
#'
#' @rdname sampler
setMethod("sampler", signature = "bdm", definition = function(object, data = list(), run = character(), init = "random", chains, iter, warmup, thin, ...) {
# initial assignments
object@data <- data
object@run <- run
# non-default sampling dimensions
if (!missing(iter)) {
object@iter <- iter
if (missing(warmup))
warmup <- floor(iter/2)
}
if (!missing(chains)) object@chains <- chains
if (!missing(thin)) object@thin <- thin
if (!missing(warmup)) {
if (warmup >= object@iter)
stop('warmup must be < iter\n')
object@warmup <- warmup
}
# check for default model
# parameterisation
pars <- getparams(object)
pars <- setdiff(pars, names(object@data))
default_pars <- setequal(pars, c('r', 'logK', 'x'))
# initial values
if (is.character(init) & default_pars) {
init.r <- getr(object)
init.logK <- getlogK(object)
init.x <- getx(object)
if (init == "fixed") {
init <- function()
list(r = init.r[['E[r]']], logK = init.logK[['E[logK]']], x = init.x[['E[x]']], sigmap = 0.05)
} else if (init == "random") {
init <- function()
list(r = rlnorm(1, init.r[['E[log(r)]']], init.r[['SD[log(r)]']]), logK = runif(1, init.logK[['min[logK]']], init.logK[['max[logK]']]), x = rbeta(length(init.x[['E[x]']]), 2, 2), sigmap = 0.05)
}
} else {
if (is.list(init)) {
init <- function() init
message("Intial values supplied on input")
}
}
# number of posterior samples
object@nsamples <- ((object@iter - object@warmup) * object@chains)/object@thin
# mcmc-sampling using rstan
stanfit_object <- suppressWarnings(sampling(object, data = object@data, init = init, iter = object@iter, chains = object@chains, warmup = object@warmup, thin = object@thin, ...))
# record initial values for each chain
object@inits <- stanfit_object@inits
# extract traces
object@trace <- extract(stanfit_object)
object@trace_array <- extract(stanfit_object,permuted = FALSE,inc_warmup = TRUE)
# check
if (dim(object@trace[['r']]) != object@nsamples) {
object@nsamples <- dim(object@trace[['r']])
}
# return
return(object)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.