R/hpd.R In BayesX: R Utilities Accompanying the Software Package BayesX

Documented in hpdhpd.coda

hpd <- function(data, alpha = 0.05, ...) {
if (!is.data.frame(data))
nc <- NCOL(data)
data <- matrix(unlist(data), ncol = nc)
data <- , -1]
hpd.target <- function(boundary, alpha, fdens, ...) {
# Constraint 1: log f(lower|x) = log f(upper|x) -- operate on log scale (more
# stable)
f1 <- log(fdens(boundary[1], ...)) - log(fdens(boundary[2], ...))

# Contraint 2: P(param \in C) = 1-alpha (use integrate -- more stable than
# ppost) f2 <- ppost(boundary[1],...) - ppost(boundary[2]) Important: the dpost
# function has to be vectorized
f2 <- integrate(fdens, lower = boundary[1], upper = boundary[2], ...)$val - (1 - alpha) # Convert multivariate root finding to minimization problem (this can be # numerically instable) return(f1^2 + f2^2) } dens <- apply(data, FUN = density, MARGIN = 2) fdens <- lapply(dens, FUN = function(obj) { splinefun(obj$x, obj$y) }) start <- apply(data, MARGIN = 2, FUN = quantile, c(0.5 * alpha, 1 - 0.5 * alpha)) opt <- function(j, ...) { optim(, j], hpd.target, alpha = alpha, fdens = fdens[[j]], ...)$par
}
res <- lapply(1:(nc - 1), FUN = opt)
res <- matrix(unlist(res), ncol = 2, byrow = TRUE)
res <- data.frame(lower = res[, 1], upper = res[, 2])
return(res)
}

hpd.coda <- function(data, alpha = 0.05) {
if (!is.data.frame(data))
nc <- NCOL(data)
data <- matrix(unlist(data), ncol = nc)
data <- , -1]
prob <- 1 - alpha
data <- as.mcmc(data)
res <- HPDinterval(data, prob = prob)

return(data.frame(lower = res[, 1], upper = res[, 2]))
}

Try the BayesX package in your browser

Any scripts or data that you put into this service are public.

BayesX documentation built on Aug. 24, 2019, 9:03 a.m.