R/summary.R

In bayesforecast: Bayesian Time Series Modeling with Stan

#' Summary method for a varstan object
#'
#' Summaries of parameter estimates and MCMC convergence diagnostics
#' (Monte Carlo error, effective sample size, Rhat).
#'
#' @param object A varstan object.
#' @param robust A boolean value, if its \code{TRUE} it returns the median of the posterior distribution,
#' And if its \code{FALSE} it returns the mean, by default is the \code{FALSE} value
#' @param prob A number \eqn{p \in (0,1)}{p (0 < p < 1)} indicating the desired
#'   probability mass to include in the intervals. The default is to report
#'   \eqn{90}% intervals (\code{prob=0.9}) rather than the traditionally used
#'   \eqn{95}%.
#' @param ... Further arguments passed to  \code{summary}.
#'
#' @author Asael Alonzo Matamoros.
#'
#' @return
#' A data.frame with the posterior mean, standard error, credible intervals, effective sample
#' size (ess),and Rhat for all the model parameters in a varstan model, if \code{robust} is \code{TRUE}
#' then the posterior mean and standard error, are replaced by the posterior mean and MAD.
#'
#' @import rstan
#'
#' @method summary varstan
#' @export
#'
summary.varstan = function(object,robust = FALSE,prob = 0.90,...){

  if(!is.varstan(object))
    stop("The current object is not a varstan class")

  conf <- 1- ((1 - prob) / 2)
  qq = c(1-conf,conf)
  par1 = get_parameters(object)
  post = as.data.frame(rstan::extract(object$stanfit,par1, permuted = TRUE) )
  sum = as.data.frame(t(apply(post,2,my_sum,robust,conf)))

  if(robust)
    names(sum) = c("median","mad", paste0(qq[1]*100,"%"), paste0(qq[2]*100,"%"),"ess","Rhat")
  else
    names(sum) = c("mean","se", paste0(qq[1]*100,"%"), paste0(qq[2]*100,"%"),"ess","Rhat" )

  return(sum)
}