#' Plot Posterior Distributions
#'
#' This function plot posterior distributions of the parameters.
#' @param chains Gibbs sampling chains
#' @param trueValues numeric vector of true values
#' @param outline logical value whether showing outliers
#' @keywords Plot Posterior Distributions
#' @export
#' @import graphics
#' @import stats
#' @examples
#'
#' ssSim <- phenoSim(nSites = 2, #number of sites
#' nTSet = 30, #number of Time steps
#' beta = c(1, 2), #beta coefficients
#' sig = .01, #process error
#' tau = .1, #observation error
#' plotFlag = TRUE, #whether plot the data or not
#' miss = 0.05, #fraction of missing data
#' ymax = c(6, 3) #maximum of saturation trajectory
#' )
#'
#' ssOut <- fitCDM(x = ssSim$x, #predictors
#' nGibbs = 200,
#' nBurnin = 100,
#' z = ssSim$z,#response
#' connect = ssSim$connect, #connectivity of time data
#' quiet=TRUE)
#'
#' summ <- getGibbsSummary(ssOut, burnin = 100, sigmaPerSeason = FALSE)
#'
#' colMeans(summ$ymax)
#' colMeans(summ$betas)
#' colMeans(summ$tau)
#' colMeans(summ$sigma)
#'
#' par(mfrow = c(1,3), oma = c(1,1,3,1), mar=c(2,2,0,1), font.axis=2)
#'
#' plotPost(chains = ssOut$chains[,c("beta.1", "beta.2")], trueValues = ssSim$beta)
#' plotPost(chains = ssOut$chains[,c("ymax.1", "ymax.2")], trueValues = ssSim$ymax)
#' plotPost(chains = ssOut$chains[,c("sigma", "tau")], trueValues = c(ssSim$sig, ssSim$tau))
#'
#' mtext('Posterior distributions of the parameters', side = 3, outer = TRUE, line = 1, font = 2)
#' legend('topleft', legend = c('posterior', 'true value'), col = c('black', 'red'),
#' lty = 1, bty = 'n', cex=1.5, lwd =2)
#'
plotPost <- function(chains, trueValues = NULL, outline = FALSE){
myboxplot(chains, outline = outline, ylim = range(range(chains, na.rm = TRUE), trueValues, na.rm = TRUE), statsParam = c(.005,.25,.5,.75,.995))
if(!is.null(trueValues)){
xs <- 1:length(trueValues)
segments(x0 = xs - 0.35, y0 = trueValues, x1 = xs + 0.35, y1 = trueValues, col = 'red', lwd=2)
}
}
#This function replaces boxplot.stats --- it uses the middle 50% and
#middle 90% instead of middle 50%(IQR) and 1.5*IQR.
myboxplot.stats <- function (x, coef = NULL,
do.conf = TRUE, do.out =TRUE,
statsParam = c(.005,.25,.5,.75,.995))
{
if(length(statsParam)!=5)
{
warning('statsParams was replaced by c(.025,.25,.5,.75,.975)!')
statsParam <- c(.025,.25,.5,.75,.975)
}
nna <- !is.na(x)
n <- sum(nna)
# stats <- quantile(x, c(.05,.25,.5,.75,.95), na.rm = TRUE)
# stats <- quantile(x, c(.025,.25,.5,.75,.975), na.rm = TRUE)
stats <- quantile(x, statsParam, na.rm = TRUE)
# iqr <- diff(stats[c(2, 4)])
out <- x < stats[1] | x > stats[5]
conf <- if (do.conf)
stats[3] + c(-1.58, 1.58) * diff(stats[c(2, 4)])/sqrt(n)
list(stats = stats, n = n, conf = conf, out = x[out & nna])
}
myboxplot <- function (x, ..., range = 1.5, width = NULL, varwidth = FALSE,
statsParam = c(.025,.25,.5,.75,.975),
notch = FALSE, outline = TRUE, names, plot = TRUE, border = par("fg"),
col = NULL, log = "", pars = list(boxwex = 0.8, staplewex = 0.5,
outwex = 0.5), horizontal = FALSE, add = FALSE, at = NULL)
{
args <- list(x, ...)
namedargs <- if (!is.null(attributes(args)$names))
attributes(args)$names != ""
else rep_len(FALSE, length(args))
groups <- if (is.list(x))
x
else args[!namedargs]
if (0L == (n <- length(groups)))
stop("invalid first argument")
if (length(class(groups)))
groups <- unclass(groups)
if (!missing(names))
attr(groups, "names") <- names
else {
if (is.null(attr(groups, "names")))
attr(groups, "names") <- 1L:n
names <- attr(groups, "names")
}
cls <- sapply(groups, function(x) class(x)[1L])
cl <- if (all(cls == cls[1L]))
cls[1L]
else NULL
for (i in 1L:n) groups[i] <- list(myboxplot.stats(unclass(groups[[i]]),
range,statsParam = statsParam))
stats <- matrix(0, nrow = 5L, ncol = n)
conf <- matrix(0, nrow = 2L, ncol = n)
ng <- out <- group <- numeric(0L)
ct <- 1
for (i in groups) {
stats[, ct] <- i$stats
conf[, ct] <- i$conf
ng <- c(ng, i$n)
if ((lo <- length(i$out))) {
out <- c(out, i$out)
group <- c(group, rep.int(ct, lo))
}
ct <- ct + 1
}
if (length(cl) && cl != "numeric")
oldClass(stats) <- cl
z <- list(stats = stats, n = ng, conf = conf, out = out,
group = group, names = names)
if (plot) {
if (is.null(pars$boxfill) && is.null(args$boxfill))
pars$boxfill <- col
do.call("bxp", c(list(z, notch = notch, width = width,
varwidth = varwidth, log = log, border = border,
pars = pars, outline = outline, horizontal = horizontal,
add = add, at = at), args[namedargs]))
invisible(z)
}
else z
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.