summarise.samples <- function(samples, columns=NULL, quantiles=0.5, exceedences=NULL)
{
#### Compute exceedence probabilities and posterior quantiles
#### for a set of given samples
#### Check for initial errors
if(class(samples) != "mcmc") stop("The samples object is not a mcmc type.", call.=FALSE)
if(!is.null(exceedences))
{
if(!is.numeric(exceedences) ) stop("the exceedences argument is not numeric.", call.=FALSE)
}else
{
}
if(!is.null(quantiles))
{
if(!is.numeric(quantiles)) stop("the quantiles argument is not numeric.", call.=FALSE)
if(min(quantiles) < 0) stop("some of the quantiles are less than zero.", call.=FALSE)
if(max(quantiles) > 1) stop("some of the quantiles are greater than one.", call.=FALSE)
}else
{
}
if(is.null(columns))
{
N.all <- ncol(samples)
columns <- 1:N.all
}else
{
N.total <- ncol(samples)
N.all <- length(columns)
if(!is.numeric(columns)) stop("the columns argument is not numeric.", call.=FALSE)
if(sum(floor(columns)==ceiling(columns)) < N.all) stop("the columns argument is not all integers.", call.=FALSE)
if(min(columns)<1) stop("some of the columns chosen are less than 1.", call.=FALSE)
if(max(columns)>N.total) stop("some of the columns chosen are greater than the size of the samples object.", call.=FALSE)
}
#### Compute the summaries
n.sample <- nrow(samples)
if(missing(exceedences))
{
N.exceedence <- 0
}else
{
N.exceedence <- length(exceedences)
}
if(missing(quantiles))
{
N.quantiles <- 0
}else
{
N.quantiles <- length(quantiles)
}
if(N.exceedence>0 & N.quantiles>0)
{
posterior.exceedence <- array(NA, c(N.all,N.exceedence))
posterior.quantile <- array(NA, c(N.all,N.quantiles))
colnames(posterior.exceedence) <- exceedences
colnames(posterior.quantile) <- quantiles
for(j in 1:N.all)
{
mcmc <- samples[ ,columns[j]]
posterior.quantile[j, ] <- quantile(mcmc, quantiles)
for(k in 1:N.exceedence)
{
posterior.exceedence[j, k] <- length(which(mcmc>exceedences[k]))/n.sample
}
}
}else if(N.exceedence==0 & N.quantiles>0)
{
posterior.exceedence <- NULL
posterior.quantile <- array(NA, c(N.all,N.quantiles))
colnames(posterior.quantile) <- quantiles
for(j in 1:N.all)
{
mcmc <- samples[ ,columns[j]]
posterior.quantile[j, ] <- quantile(mcmc, quantiles)
}
}else if(N.exceedence>0 & N.quantiles==0)
{
posterior.exceedence <- array(NA, c(N.all,N.exceedence))
posterior.quantile <- NULL
colnames(posterior.exceedence) <- exceedences
for(j in 1:N.all)
{
mcmc <- samples[ ,columns[j]]
for(k in 1:N.exceedence)
{
posterior.exceedence[j, k] <- length(which(mcmc>exceedences[k]))/n.sample
}
}
}else
{
posterior.quantile <- NULL
posterior.exceedence <- NULL
}
results <- list(quantiles=posterior.quantile, exceedences=posterior.exceedence)
return(results)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.