R/hpdparameter.R

In QAEnsemble: Ensemble Quadratic and Affine Invariant Markov Chain Monte Carlo

#' Highest Posterior Density (HPD) for a parameter
#'
#'@description This function returns the upper and lower bound of the Highest
#'Posterior Density (HPD) for a parameter based on the Chen-Shao Highest
#'Posterior Density (HPD) Estimation Algorithm found in the book by Chen, Shao,
#'and Ibrahim (2000). (The smallest 95% credible interval will be given by the
#'HPD using alpha = 0.05)
#'
#' @param parameter_MCMC a vector of the parameter samples for a single
#' estimated parameter
#' @param alpha 100(1 - alpha)% credible interval with the default value as
#' alpha = 0.05
#'
#' @return A vector is returned that contains the lower and upper bound of the
#' Highest Posterior Density (HPD) for a parameter (this will be the smallest
#' 95% credible interval using alpha = 0.05)
#' @export
#'
#' @references Chen M, Shao Q, Ibrahim JG (2000) Monte Carlo Methods in Bayesian
#' Computation. New York-New York: Springer-Verlag.
#'
#' @examples x_parameter = rnorm(75, mean = 0, sd = 1)
#'
#' hpdparameter(x_parameter, 0.05)
hpdparameter <- function(parameter_MCMC, alpha = 0.05)
{
  parameter_MCMC_ordered = sort(parameter_MCMC)

  n = length(parameter_MCMC)

  R_j = matrix(NA, nrow = 2, ncol = (n - floor((1 - alpha)*n)))

  for(j in 1:(n - floor((1 - alpha)*n)))
  {
    R_j[1,j] = parameter_MCMC_ordered[j]
    R_j[2,j] = parameter_MCMC_ordered[j + floor((1 - alpha)*n)]
  }

  R_j_diff = matrix(NA, nrow = 1, ncol = (n - floor((1 - alpha)*n)))

  for(j in 1:(n - floor((1 - alpha)*n)))
  {
    R_j_diff[1, j] = R_j[2,j] - R_j[1,j]
  }

  index_use = which.min(R_j_diff)

  return(c(R_j[1,index_use], R_j[2, index_use]))

}