R/fit.R

In gppm: Gaussian Process Panel Modeling

#' Generic Method For Fitting a model
#'
#' Generic method for fitting a model.
#' @param gpModel a model.
#' @param ... additional arguments.
#' @return A fitted model
#' @seealso \code{\link{fit.GPPM}}
#' @export
fit <-  function(gpModel,...) {
  UseMethod("fit")
}

#' Fit a Gaussian process panel model
#'
#' This function is used to fit a Gaussian process panel model,
#' which has been specified fit using \code{\link{gppm}}.
#'
#' @param gpModel object of class GPPM. The Gaussian process panel model to be fitted.
#'
#' @param init string or named numeric vector. Used to specify the starting values for the parameters. Can either be the string 'random' (default) or a numeric vector startVal of starting values. Which value belongs to which parameter is determined by the names attribute of startVal. See also the example.
#'
#' @param useOptimizer boolean. Should the optimizer be used or not? For false the (possibly random) starting values are returned as the maximum likelihood estimates.
#'
#' @param verbose boolean. Print diagnostic output?
#'
#' @param hessian boolean. Compute the hessian at the maximum likelihood estimate?
#'
#' @param ... additional arguments (currently not used).
#'
#' @return A fitted Gaussian process panel model, which is an object of class 'GPPM'.
#' @examples
#' \donttest{
#' #regular usage
#' data("demoLGCM")
#' lgcm <- gppm('muI+muS*t','varI+covIS*(t+t#)+varS*t*t#+(t==t#)*sigma',
#'         demoLGCM,'ID','y')
#' lgcmFit <- fit(lgcm)
#'
#' #starting values as ML results
#' startVals <- c(10,1,10,3,10,1)
#' names(startVals) <- pars(lgcm)
#' lgcmFakeFit <- fit(lgcm,init=startVals,useOptimizer=FALSE)
#' stopifnot(identical(startVals,coef(lgcmFakeFit)))
#' }
#' @export
fit.GPPM <-  function(gpModel,init='random',useOptimizer=TRUE,verbose=FALSE,hessian=TRUE,...) {
  if (useOptimizer){
    iter<- 10000
    algorithm <- 'LBFGS' #default
  }
  else{
    iter <- 0
    algorithm <- 'Newton'#because for some reason LBFSGS still does an iteration even if iter==0
  }
  if (!init[1]=='random'){
    validate_simulate(gpModel,init)
    init = as.list(init)
  }
  if(verbose){
    outf <- eval
  }else{
    outf <- utils::capture.output
  }
  outf(gpModel$stanOut <- rstan::optimizing(gpModel$stanModel,gpModel$dataForStan,hessian = hessian,iter=iter,init=init,algorithm=algorithm,as_vector=FALSE))
  gpModel$fitRes <- extractFitRes(gpModel$stanOut,gpModel$parsedModel,gpModel$dataForStan[c('nPer','nTime','maxTime','Y')])
  gpModel
}