R/sampler.R

In cttedwards/bdm: Bayesian biomass dynamics model

#'
#' 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)
})