#' This function is used to predict the age of new stratigraphic points using the output from an ageModel run.
#'
#' @param model Output of the \code{ageModel} function
#' @param newPositions Vector of new stratigraphic positions to predict the age of. New positions must be within the bounds of the orignal model run
#' @param ids optional vector of names for new positions. Defaults to 1, 2, 3, ...
#' @param newPositionThicknesses Vector of stratigraphic uncertanties for each position. Specified as half thicknesses. Must be the same length and given in the same order as \code{newPositions}. Default is no positional uncertanty
#' @return HDI = 95 percent Highest Density Interval for each \code{newPosition}
#' @return raw = Age predictions for each \code{newPosition} for each MCMC run.
#' @export
agePredict <- function(model,
newPositions,
ids = 1:length(newPositions),
newPositionThicknesses = rep(0, length(newPositions)),
probability = 0.95){
## this function is used to predict the age of new points using a modified bchron output model
## INPUTS
## model = output of the modified bchron age-depth model
## newPositions = new stratigraphic points to predict the age of.
## The new points must fall within the original stratigraphic range of the model
## OUTPUTS
## ConfInt = 95% Confidence interval for each point
## raw = age predictions for each point for each iteration in the MCMC chain
pb <- utils::txtProgressBar(min = 1,
max = ncol(model$model),
style = 3,
width = 40,
char = '<>') # create a progress bar
x <- model$predictPositions
n <- length(newPositions)
predictStore <- matrix(nrow = ncol(model$model),
ncol = n)
for(i in 1:ncol(model$model)){
currPositions <- runif(n,
newPositions - newPositionThicknesses,
newPositions + newPositionThicknesses)
utils::setTxtProgressBar(pb, i) # set the progress bar
y <- model$model[, i]
f <- approxfun(x, y)
predictStore[i, ] <- t(f(currPositions))
}
# remove burn burn-in
predictStore <- predictStore[model$burn: ncol(model$model), ]
predictStore <- as.data.frame(predictStore)
HDI <- t(apply(predictStore, 2, quantile, c((1 - probability) / 2, 0.5, (1 + probability) / 2), na.rm = T))
HDI <- cbind(newPositions, HDI)
HDI <- data.frame(HDI)
HDI$ids <- ids
colnames(HDI) <- c('newPositions',
as.character(c((1 - probability) / 2, 0.5, (1 + probability) / 2)),
'ids')
predictStore <- data.frame(predictStore)
return(list(HDI = HDI,
raw = predictStore))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.