#' Runs the Compound Poisson-Gamma chronology model of Haslett and Parnell (2008)
#'
#' Fits a non-parametric chronology model to age/position data according to the Compound Poisson-Gamma model defined by Haslett and Parnell (2008) <DOI:10.1111/j.1467-9876.2008.00623.x>. This version uses a slightly modified Markov chain Monte Carlo fitting algorithm which aims to converge quicker and requires fewer iterations. It also a slightly modified procedure for identifying outliers
#'
#' @inheritParams BchronCalibrate
#'
#' @param positionThicknesses Thickness values for each of the positions. The thickness value should be the full thickness value of the slice. By default set to zero.
#' @param outlierProbs A vector of prior outlier probabilities, one for each age. Defaults to 0.01
#' @param predictPositions A vector of positions (e.g. depths) at which predicted age values are required. Defaults to a sequence of length 100 from the top position to the bottom position
#' @param artificialThickness Amount to add to the thickness values in the case of equal positions with no \code{positionThicknesses}. Bchron may fail if \code{positionThicknesses} are zero and some positions are repeated. This value is added on to the zero thicknesses (only in the case of repeated positions) to stop this failure.
#' @param iterations The number of iterations to run the procedure for
#' @param burn The number of starting iterations to discard
#' @param thin The step size for every iteration to keep beyond the burn-in
#' @param extractDate The top age of the core. Used for extrapolation purposes so that no extrapolated ages go beyond the top age of the core. Defaults to the current year
#' @param maxExtrap The maximum number of extrapolations to perform before giving up and setting the predicted ages to NA. Useful for when large amounts of extrapolation are required, i.e. some of the \code{predictPositions} are a long way from the dated positions
#' @param thetaStart A set of starting values for the calendar ages estimated by Bchron. If NULL uses a function to estimate the ages. These should be in the same units as the posterior ages required. See example below for usage.
#' @param thetaMhSd The Metropolis-Hastings standard deviation for the age parameters
#' @param muMhSd The Metropolis-Hastings standard deviation for the Compound Poisson-Gamma mean
#' @param psiMhSd The Metropolis-Hastings standard deviation for the Compound Poisson-Gamma scale
#' @param ageScaleVal A scale value for the ages. \code{Bchronology} works best when the ages are scaled to be approximately between 0 and 100. The default value is thus 1000 for ages given in years.
#' @param positionEps A small value used to check whether simulated positions are far enough apart to avoid numerical underflow errors. If errors occur in model runs (e.g. \code{missing value where TRUE/FALSE needed} increase this value)
#' @param positionNormalise Whether to normalise the position values. \code{Bchronology} works best when the positions are normalised to be between 0 and 1 The default value is \code{TRUE}
#'
#' @details
#' The \code{Bchronology} function fits a compound Poisson-Gamma distribution to the increments between the dated levels. This involves a stochastic linear interpolation step where the age gaps are Gamma distributed, and the position gaps are Exponential. Radiocarbon and non-radiocarbon dates (including outliers) are updated within the function also by MCMC.
#'
#' @useDynLib Bchron
#'
#' @return A list of class \code{BchronologyRun} which include elements:
#' \item{theta}{The posterior estimated values of the ages}
#' \item{phi}{The posterior estimated outlier values (1=outlier, 2=not outlier). The means of this parameter give the posterior estimated outlier probabilities}
#' \item{mu}{The posterior values of the Compound Poisson-Gamma mean}
#' \item{psi}{The posterior values of the Compound Poisson-Gamma scale}
#' \item{thetaPredict}{The posterior estimated ages for each of the values in predictPosition}
#' \item{predictPositions}{The positions at which estimated ages were required}
#' \item{calAges}{The calibrated ages as output from \code{\link{BchronCalibrate}}}
#' \item{inputVals}{All of the input values to the \code{Bchronology} run}
#'
#' @references
#' Haslett, J., and Parnell, A. C. (2008). A simple monotone process with application to radiocarbon-dated depth chronologies. Journal of the Royal Statistical Society, Series C, 57, 399-418. DOI:10.1111/j.1467-9876.2008.00623.x
#' Parnell, A. C., Haslett, J., Allen, J. R. M., Buck, C. E., and Huntley, B. (2008). A flexible approach to assessing synchroneity of past events using Bayesian reconstructions of sedimentation history. Quaternary Science Reviews, 27(19-20), 1872-1885. DOI:10.1016/j.quascirev.2008.07.009
#'
#' @seealso
#' \code{\link{BchronCalibrate}}, \code{\link{BchronRSL}}, \code{\link{BchronDensity}}, \code{\link{BchronDensityFast}}
#' @examples
#' \donttest{
#' # Data from Glendalough
#' data(Glendalough)
#'
#' # Run in Bchronology - all but first age uses intcal20
#' GlenOut <- with(
#' Glendalough,
#' Bchronology(
#' ages = ages,
#' ageSds = ageSds,
#' calCurves = calCurves,
#' positions = position,
#' positionThicknesses = thickness,
#' ids = id,
#' predictPositions = seq(0, 1500, by = 10)
#' )
#' )
#'
#' # Summarise it a few different ways
#' summary(GlenOut) # Default is for quantiles of ages at predictPosition values
#' summary(GlenOut, type = "convergence") # Check model convergence
#' summary(GlenOut, type = "outliers") # Look at outlier probabilities
#'
#' # Predict for some new positions
#' predictAges <- predict(GlenOut,
#' newPositions = c(150, 725, 1500),
#' newPositionThicknesses = c(5, 0, 20)
#' )
#'
#' # Plot the output
#' plot(GlenOut) +
#' ggplot2::labs(
#' title = "Glendalough",
#' xlab = "Age (cal years BP)",
#' ylab = "Depth (cm)"
#' )
#'
#' # If you need to specify your own starting values
#' startingAges <- c(0, 2000, 10000, 11000, 13000, 13500)
#' GlenOut <- with(
#' Glendalough,
#' Bchronology(
#' ages = ages,
#' ageSds = ageSds,
#' calCurves = calCurves,
#' positions = position,
#' positionThicknesses = thickness,
#' ids = id,
#' predictPositions = seq(0, 1500, by = 10),
#' thetaStart = startingAges
#' )
#' )
#' }
#'
#' @export
Bchronology <- function(ages,
ageSds,
positions,
positionThicknesses = rep(
0,
length(ages)
),
calCurves = rep(
"intcal20",
length(ages)
),
ids = NULL,
outlierProbs = rep(0.01, length(ages)),
predictPositions = seq(min(positions),
max(positions),
length = 100
),
pathToCalCurves = system.file("data", package = "Bchron"),
artificialThickness = 0.01,
allowOutside = FALSE,
iterations = 10000,
burn = 2000,
thin = 8,
extractDate = 1950 - as.numeric(format(Sys.time(), "%Y")),
maxExtrap = 1000,
thetaStart = NULL,
thetaMhSd = 0.5,
muMhSd = 0.1,
psiMhSd = 0.1,
ageScaleVal = 1000,
positionEps = 1e-5,
positionNormalise = TRUE) {
# Notation:
# theta are the calibrated ages of ages 1 to n (not necessarily radiocarbon)
# phi are the outlier indicators (1=TRUE or 0=FALSE) for date i
# mu,psi are the Compound Poisson-Gamma parameters controlling sedimentation
# Run data checks ---------------------------------------------------------
# Run the Bchron check function in case things have gone wrong
BchronCheck(
ages = ages,
ageSds = ageSds,
positions = positions,
positionThicknesses = positionThicknesses,
calCurves = calCurves,
ids = ids,
outlierProbs = outlierProbs,
predictPositions = predictPositions,
pathToCalCurves = pathToCalCurves,
artificialThickness = artificialThickness,
allowOutside = allowOutside,
iterations = iterations,
thetaStart = thetaStart,
burn = burn,
thin = thin,
extractDate = extractDate,
maxExtrap = maxExtrap,
thetaMhSd = thetaMhSd,
muMhSd = muMhSd,
psiMhSd = psiMhSd,
ageScaleVal = ageScaleVal,
positionEps = positionEps,
positionNormalise = positionNormalise,
type = "Bchronology"
)
# Re-normalise positions --------------------------------------------------
originalPositions <- positions
originalPositionThicknesses <- positionThicknesses
originalPredictPositions <- predictPositions
if (positionNormalise) {
positionRange <- diff(range(positions))
positions <- (positions - min(originalPositions)) / positionRange
positionThicknesses <- positionThicknesses / positionRange
predictPositionsRescaled <- (predictPositions - min(originalPositions)) / positionRange
} else {
predictPositionsRescaled <- predictPositions
}
# Check positions don't overlap
n <- length(ages)
depthLow <- positions - 0.5 * positionThicknesses
depthHigh <- positions + 0.5 * positionThicknesses
# Function to create decent starting positions
getCurrPositions <- function(pos, posThick, posEps) {
n <- length(pos)
badPositions <- TRUE
while (badPositions) {
currPositions <- stats::runif(
n,
pos - 0.5 * posThick,
pos + 0.5 * posThick
)
do <- order(currPositions)
diffPosition <- diff(currPositions[do])
if (all(diffPosition > posEps)) badPositions <- FALSE
}
return(currPositions)
}
# Check order of positions ------------------------------------------------
# Check positions are in order
o <- order(positions)
if (any(positions[o] != positions)) {
warning("positions not given in order - re-ordering")
ages <- ages[o]
ageSds <- ageSds[o]
positions <- positions[o]
positionThicknesses <- positionThicknesses[o]
calCurves <- calCurves[o]
ids <- ids[o]
outlierProbs <- outlierProbs[o]
}
# Calibrate dates ---------------------------------------------------------
# First thing to do is to calibrate the ages with the two different degrees of freedom
dfs <- c(1, 100) # corresponding to phi equals 1 and 0 respectively
x.df1 <- BchronCalibrate(
ages = ages,
ageSds = ageSds,
calCurves = calCurves,
positions = originalPositions,
ids = ids,
pathToCalCurves = pathToCalCurves,
allowOutside = allowOutside,
eps = 0,
dfs = rep(dfs[1], length(ages))
)
x.df2 <- BchronCalibrate(
ages = ages,
ageSds = ageSds,
calCurves = calCurves,
positions = originalPositions,
ids = ids,
pathToCalCurves = pathToCalCurves,
allowOutside = allowOutside,
eps = 0,
dfs = rep(dfs[2], length(ages))
)
# Finally - need to make sure that all calibrated dates have the same ageGrids
# This is to avoid the situation where the depth order swaps and you end up
# with ages that run outside their calibration curve ranges
# The densities outside this range are set to zero
range1 <- range(lapply(x.df1, "[", "ageGrid"))
range2 <- range(lapply(x.df2, "[", "ageGrid"))
masterAgeGrid1 <- seq(min(range1), max(range1), by = 1)
masterAgeGrid2 <- seq(min(range2), max(range2), by = 1)
for (j in 1:length(x.df1)) {
if (!setequal(x.df1[[j]]$ageGrid, masterAgeGrid1)) {
currAgeGrid <- x.df1[[j]]$ageGrid
currDensities <- x.df1[[j]]$densities
x.df1[[j]]$densities <- rep(0, length(masterAgeGrid1))
x.df1[[j]]$densities[match(x.df1[[j]]$ageGrid, masterAgeGrid1)] <- currDensities
x.df1[[j]]$ageGrid <- masterAgeGrid1
}
if (!setequal(x.df2[[j]]$ageGrid, masterAgeGrid2)) {
currAgeGrid <- x.df2[[j]]$ageGrid
currDensities <- x.df2[[j]]$densities
x.df2[[j]]$densities <- rep(0, length(masterAgeGrid2))
x.df2[[j]]$densities[match(x.df2[[j]]$ageGrid, masterAgeGrid2)] <- currDensities
x.df2[[j]]$ageGrid <- masterAgeGrid2
}
}
# Get starting positions --------------------------------------------------
# If it has position thicknesses for all then use these to create some current positions
if (any(positionThicknesses == 0)) {
if (length(positionThicknesses[duplicated(positions)]) > 0) {
message("Some positionThicknesses are zero for identical positions. artificialThickness has been added so that the model can attempt to run. If the model still fails then increase the value of artificialThickness further.")
positionThicknesses[duplicated(positions)] <- positionThicknesses[duplicated(positions)] + artificialThickness
currPositions <- getCurrPositions(
positions,
positionThicknesses, positionEps
)
} else {
currPositions <- getCurrPositions(
positions,
positionThicknesses, positionEps
)
}
} else {
currPositions <- getCurrPositions(
positions,
positionThicknesses, positionEps
)
}
do <- order(currPositions)
diffPosition <- diff(currPositions[do])
# For all dates we need an offset to enable fast lookup
# The offset will be the same for all if we're using a master ageGrid
offset <- rep(-min(x.df1[[1]]$ageGrid), length = n)
# Useful C functions ------------------------------------------------------
# C function for truncated random walk
truncatedWalk <- function(old, sd, low, high) {
if (isTRUE(all.equal(low, high, tolerance = 1e-4))) {
return(list(new = low, rat = 1))
}
if (high < low) warning("truncatedWalk has lower limit higher than upper limit")
new <- .C(
"truncatedWalk",
as.double(old),
as.double(sd),
as.double(low),
as.double(high),
as.double(0)
)[5][[1]]
rat <- .C(
"truncatedRat",
as.double(old),
as.double(sd),
as.double(low),
as.double(high),
as.double(new),
as.double(0)
)[6][[1]]
if (is.nan(rat)) {
rat <- 1
} # Useful for when proposed value and new value are identical
return(list(new = new, rat = rat))
}
# C functions for tweedie
dtweediep1 <- Vectorize(function(y, p, mu, phi) {
if (phi < 0) stop("Bad phi parameter value")
if (mu < 0) stop("Bad mu parameter value")
return(.C(
"dtweediep1",
as.double(y),
as.double(p),
as.double(mu),
as.double(phi),
as.double(0)
)[5][[1]])
})
# Some C functions to do prediction and interpolation
predictInterp <- function(alpha,
lambda,
beta,
predictPositions,
diffPositionj,
currPositionsj,
currPositionsjp1,
thetaj,
thetajp1) {
return(
.C(
"predictInterp",
as.double(alpha),
as.double(lambda),
as.double(beta),
as.double(predictPositions),
as.integer(length(predictPositions)),
as.double(diffPositionj),
as.double(currPositionsj),
as.double(currPositionsjp1),
as.double(thetaj),
as.double(thetajp1),
as.double(rep(NA_real_, length(
predictPositions
))),
NAOK = TRUE
)[11][[1]]
)
}
predictExtrapUp <- function(alpha,
lambda,
beta,
predictPositions,
currPositions1,
theta1,
maxExtrap,
extractDate) {
return(
.C(
"predictExtrapUp",
as.double(alpha),
as.double(lambda),
as.double(beta),
as.double(predictPositions),
as.integer(length(predictPositions)),
as.double(currPositions1),
as.double(theta1),
as.integer(maxExtrap),
as.double(extractDate),
as.double(rep(NA_real_, length(
predictPositions
))),
NAOK = TRUE
)[10][[1]]
)
}
predictExtrapDown <- function(alpha,
lambda,
beta,
predictPositions,
currPositionsn,
thetan,
maxExtrap) {
return(
.C(
"predictExtrapDown",
as.double(alpha),
as.double(lambda),
as.double(beta),
as.double(predictPositions),
as.integer(length(predictPositions)),
as.double(currPositionsn),
as.double(thetan),
as.integer(maxExtrap),
as.double(rep(0, length(
predictPositions
))),
NAOK = TRUE
)[9][[1]]
)
}
# End of C functions
# Starting values ---------------------------------------------------------
# Written by Nathan McJames - calculate starting values
if (is.null(thetaStart)) {
thetaStart <- function(a, b, ageGrids, ageDensities) {
minbs <- rev(cummin(rev(b)))
maxas <- cummax(a)
if (any(minbs < maxas)) {
stop("Invalid starting values. Use the thetaStart argument to specify valid starting calendar ages")
} else {
raw_soln <- rep(NA, length(b))
for (u in 1:length(b)) {
currAgeGrid <- ageGrids[[u]]$ageGrid / ageScaleVal
currDensities <- ageDensities[[u]]$densities
# Remove previous values from ageGrid
currAgeGrid <- currAgeGrid[which(!currAgeGrid %in% raw_soln[1:(u - 1)])]
currDensities <- currDensities[which(!currAgeGrid %in% raw_soln[1:(u - 1)])]
# Now derive raw (unsorted solutions)
raw_soln[u] <- sample(
x = currAgeGrid[currAgeGrid >= maxas[u] & currAgeGrid <= minbs[u]],
size = 1,
prob = currDensities[currAgeGrid >= maxas[u] & currAgeGrid <= minbs[u]]
)
}
soln <- sort(raw_soln)
if (any(diff(soln) < 1e-4)) stop(paste0("Invalid starting values. Use the thetaStart argument to specify valid starting calendar ages. Bchron's initial guess was c(", paste(soln, collapse = ", "), ")"))
return(sort(soln))
}
}
ageGrids <- lapply(x.df2, "[", "ageGrid")
ageDensities <- lapply(x.df2, "[", "densities")
ranges <- do.call(rbind, lapply(ageGrids, range)) / ageScaleVal
theta <- thetaStart(ranges[, 1], ranges[, 2], ageGrids, ageDensities)
} else {
theta <- thetaStart / ageScaleVal
}
# Make sure that the thetas are in the same order as the current positions
theta[do] <- theta
# Other starting values
phi <- rep(0, length(theta))
p <- 1.2
mu <- abs(stats::rnorm(
1,
mean = mean(diff(theta)) / mean(diffPosition),
sd = 1
))
psi <- abs(stats::rnorm(1, 1, 1))
# Tranformed values (used for interpolation)
alpha <- (2 - p) / (p - 1)
# Storage -----------------------------------------------------------------
# Storage
remaining <- (iterations - burn) / thin
thetaStore <- phiStore <- matrix(ncol = length(theta), nrow = remaining)
muStore <- psiStore <- vector(length = remaining)
thetaPredict <- matrix(ncol = length(predictPositions), nrow = remaining)
# Main iteration loop -----------------------------------------------------
cat("Running Bchronology...\n")
pb <- utils::txtProgressBar(
min = 1,
max = iterations,
style = 3,
width = 60,
)
for (i in 1:iterations) {
utils::setTxtProgressBar(pb, i)
# Update positions --------------------------------------------------------
if (any(positionThicknesses > 0) &
i > 0.5 * burn & i %% thin == 0) {
# Get date order so I can preserve things if they change around
currPositions <- getCurrPositions(positions, positionThicknesses, positionEps)
do <- order(currPositions)
diffPosition <- diff(currPositions[do])
theta[do] <- sort(theta)
}
# Put things in storage ---------------------------------------------------
# If we're in the right place put things in storage
if (i > burn & i %% thin == 0) {
ind <- (i - burn) / thin
thetaStore[ind, ] <- theta * ageScaleVal
phiStore[ind, ] <- phi
muStore[ind] <- mu
psiStore[ind] <- psi
# Run interpolation/extrapolation stage
lambda <- (mu^(2 - p)) / (psi * (2 - p))
beta <- 1 / (psi * (p - 1) * (mu^(p - 1)))
# Interpolation and extrapolation -----------------------------------------
# First interpolation
for (j in 1:(n - 1)) {
# Find which positions we need to interpolate for
depthIndRange <- which(
predictPositionsRescaled >= currPositions[do[j]] &
predictPositionsRescaled <= currPositions[do[j + 1]]
)
if (length(depthIndRange) > 0) {
thetaPredict[ind, depthIndRange] <- round(
predictInterp(
alpha,
lambda,
beta,
predictPositionsRescaled[depthIndRange],
diffPosition[j],
currPositions[do[j]],
currPositions[do[j + 1]],
theta[do[j]],
theta[do[j + 1]]
),
3
)
}
}
# Extrapolate up to to top depth
if (any(predictPositionsRescaled < currPositions[1])) {
depthIndRange <- which(predictPositionsRescaled <= currPositions[1])
thetaPredict[ind, depthIndRange] <- round(
predictExtrapUp(
alpha,
lambda,
beta,
predictPositionsRescaled[depthIndRange],
currPositions[1],
theta[1],
maxExtrap,
extractDate / ageScaleVal
),
3
)
}
# Extrapolate below bottom depth
if (any(predictPositionsRescaled >= currPositions[n])) {
depthIndRange <- which(predictPositionsRescaled >= currPositions[n])
thetaPredict[ind, depthIndRange] <- round(
predictExtrapDown(
alpha,
lambda,
beta,
predictPositionsRescaled[depthIndRange],
currPositions[n],
theta[n],
maxExtrap
),
3
)
}
if (any(is.na(thetaPredict[ind, ]))) {
warning(
"NA values in predicted ages. Check you are not extrapolating too far away from dated levels. If you must run this core with these values of predictPositions, set maxExtrap to a larger value (e.g. 1000). "
)
}
}
# Update theta ------------------------------------------------------------
for (j in 1:n) {
lowerLimit <- ifelse(
j == 1,
ifelse(x.df1[[do[j]]]$calCurve == "normal",
extractDate / ageScaleVal,
min(x.df1[[do[j]]]$ageGrid) / ageScaleVal
),
theta[do[j - 1]] + 0.001
)
# upperLimit <- ifelse(j == n, 100000, theta[do[j + 1]] - 0.001)
upperLimit <- ifelse(j == n,
max(x.df1[[do[j]]]$ageGrid) / ageScaleVal,
theta[do[j + 1]] - 0.001
)
thetaNewAll <- truncatedWalk(
theta[do[j]],
thetaMhSd,
lowerLimit,
upperLimit
)
thetaNew <- round(thetaNewAll$new, 3)
# Calculate ratio
if (phi[do[j]] == 0) {
currDens <- x.df2[[do[j]]]$densities
} else {
currDens <- x.df1[[do[j]]]$densities
}
thetaNewMatch <- as.integer(thetaNew * ageScaleVal + offset[do[j]]) +
1
thetaNewLogDens <- max(log(currDens[thetaNewMatch]), -1000000,
na.rm =
TRUE
) # Get rid of NAs in case of moving beyond the calibration curve
priorNewLogDens <- ifelse(j == 1, 0, log(dtweediep1(
thetaNew - theta[do[j - 1]],
p,
mu * diffPosition[j - 1],
psi / (diffPosition[j - 1]^(p - 1))
))) + ifelse(j == n, 0, log(dtweediep1(
theta[do[j + 1]] - thetaNew,
p,
mu * (diffPosition[j]),
psi / (diffPosition[j])^(p - 1)
)))
thetaMatch <- as.integer(theta[do[j]] * ageScaleVal + offset[do[j]]) +
1
thetaLogDens <- max(log(currDens[thetaMatch]), -1000000)
priorLogDens <- ifelse(j == 1, 0, log(dtweediep1(
theta[do[j]] - theta[do[j - 1]],
p,
mu * (diffPosition[j - 1]),
psi / (diffPosition[j - 1])^(p - 1)
))) + ifelse(j == n, 0, log(dtweediep1(
theta[do[j + 1]] - theta[do[j]],
p,
mu * (diffPosition[j]),
psi / (diffPosition[j])^(p - 1)
)))
logRtheta <- thetaNewLogDens - thetaLogDens + priorNewLogDens - priorLogDens + log(thetaNewAll$rat)
if (stats::runif(1) < exp(logRtheta)) {
theta[do[j]] <- thetaNew
}
}
# Update other parameters -------------------------------------------------
# Update phi
for (j in 1:n) {
phiNew <- sample(0:1, 1)
if (phiNew != phi[do[j]]) {
if (phi[do[j]] == 0) {
currDens <- x.df2[[do[j]]]$densities
} else {
currDens <- x.df1[[do[j]]]$densities
}
thetaMatch <- as.integer(theta[do[j]] * ageScaleVal + offset[j]) +
1
thetaLogDens <- max(log(currDens[thetaMatch]), -1000000)
if (phiNew == 0) {
newDens <- x.df2[[do[j]]]$densities
} else {
newDens <- x.df1[[do[j]]]$densities
}
thetaMatch <- as.integer(theta[do[j]] * ageScaleVal + offset[j]) +
1
thetaNewLogDens <- max(log(newDens[thetaMatch]), -1000000)
logRphi <- thetaNewLogDens - thetaLogDens + stats::dbinom(phiNew, 1, outlierProbs[do[j]],
log =
TRUE
) - stats::dbinom(phi[do[j]], 1, outlierProbs[do[j]], log = TRUE)
if (stats::runif(1) < exp(logRphi)) {
phi[do[j]] <- phiNew
}
}
}
# Update mu
muNewAll <- truncatedWalk(mu, muMhSd, 1e-3, 1e3)
muNew <- muNewAll$new
logRmu <- sum(log(dtweediep1(
diff(theta[do]),
p,
muNew * diffPosition,
psi / (diffPosition)^(p - 1)
))) - sum(log(dtweediep1(
diff(theta[do]),
p,
mu * diffPosition,
psi / (diffPosition)^(p - 1)
))) + log(muNewAll$rat)
if (stats::runif(1) < exp(logRmu)) {
mu <- muNew
}
# Update psi
psiNewAll <- truncatedWalk(psi, psiMhSd, 1e-3, 1e3)
psiNew <- psiNewAll$new
logRpsi <- sum(log(dtweediep1(
diff(theta[do]),
p,
mu * diffPosition,
psiNew / (diffPosition)^(p - 1)
))) - sum(log(dtweediep1(
diff(theta[do]),
p,
mu * diffPosition,
psi / (diffPosition)^(p - 1)
))) + log(psiNewAll$rat)
if (stats::runif(1) < exp(logRpsi)) {
psi <- psiNew
}
}
# Collect up and return ---------------------------------------------------
# Collect up the input values to return
inputVals <- list(
ages = ages,
ageSds = ageSds,
positions = originalPositions,
positionThicknesses = originalPositionThicknesses,
positionsNormalised = positions,
positionThicknessesNormalised = positionThicknesses,
calCurves = calCurves,
ids = ids,
outlierProbs = outlierProbs,
predictPositions = predictPositions,
pathToCalCurves = pathToCalCurves,
artificialThickness = artificialThickness,
iterations = iterations,
burn = burn,
thin = thin,
extractDate = extractDate,
maxExtrap = maxExtrap,
thetaMhSd = thetaMhSd,
muMhSd = muMhSd,
psiMhSd = psiMhSd,
ageScaleVal = ageScaleVal,
positionNormalise = positionNormalise
)
# Return everything
out <- list(
theta = thetaStore,
phi = phiStore,
mu = muStore,
psi = psiStore,
thetaPredict = ageScaleVal * thetaPredict,
predictPositions = originalPredictPositions,
calAges = x.df2,
positions = originalPositions,
extractDate = extractDate,
ageScaleVal = ageScaleVal,
positionNormalise = positionNormalise,
inputVals = inputVals
)
class(out) <- "BchronologyRun"
cat("\nRun completed!\n")
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.