R/hdr.R

In twolodzko/twextras: twextras

#' Calculate Highest Density Interval from vector of probabilities
#' 
#' @param x        numeric vector of probabilities or density estimates.
#' @param prob     probability coverage.
#' @param method   method of calculating HDR to be used: \code{"quantile"} for
#'                 using quantiles as suggested by Hyndman (1996), \code{"grid"} for 
#'                 grid-based method suggested by Kruschke (2011).
#' @param simulate if \code{TRUE} bootstrap samples are drawn from \code{x} for
#'                 the \code{"quantile"} method.
#' @param nsim     number of samples to be drawn if \code{simulate=TRUE}.
#' 
#' @references 
#' Hyndman, R.J. (1996) Computing and graphing highest density regions.
#' American Statistician, 50, 120-126.
#' 
#' @references 
#' Kruschke, J.K. (2011). Doing Bayesian data analysis: A Tutorial with R and BUGS.
#' Elsevier Science/Academic Press.
#' 
#' @importFrom stats quantile
#' @export

hdr <- function(x, prob = 0.95, method = c("quantile", "grid"),
                simulate = TRUE, nsim = 5000) {
  
  method <- match.arg(method)
  
  if (method == "grid") {
    
    const <- sum(x)
    if (const != 1)
      x <- x/const
    x <- sort(x, decreasing = TRUE)
    csx <- cumsum(x)
    idx <- vapply(prob, function(p) which(csx >= p)[1], numeric(1))
    height <- x[idx]*const
    
  } else {
    
    if (simulate)
      x <- sample(x, nsim, replace = TRUE, prob = x)
    height <- quantile(x, probs = 1-prob)
    
  }
  
  names(height) <- stats:::format_perc(prob)
  return(height)
}