R/calibrate.r

In rmvegasm/schron: Bayesian chronologies for multiple sedimentary sequences.

#' @title Calibrate 14C age determinations
#'
#' @param age a numeric vector containing the mean 14C age for the samples.
#'
#' @param error a numeric vector containing the error for each 14C
#' determination.
#'
#' @param id an optional character vector of labels for each 14C age, used for
#' names in the output. 
#'
#' @param curve a character string giving the name of the calibration curve:
#' - 'intcal20'
#' - 'shcal20'
#' - 'marine20'
#' - 'intcal13'
#' - 'shcal13'
#' - 'marine13'
#' - 'intcal09'
#'
#' @param eps the smallest probability to be considered as something. Values
#' under this number will not be included in the `pmf` output.
#'
#' @details Calibration is done on a yearly defined grid. The probability for
#' each year is computed by the same method used by OxCal: if the 14C date is
#' r_m with a measurement error σ_m the variance distribution v_i and
#' the probability distribution p_i are given by:
#'
#' v_i = (r_m - r_i)^2 / (σ_m^2 + σ_i^2)
#'
#' p_i ~ exp(-v_i / 2) / sqrt(σ_m^2 + σi^2)
#'
#' A Probability Mass Function (PMF) is constructed for those values with
#' p_i >= eps.
#'
#' @return a list of `pmf` objects.
#' @export
calibrate <- function (age, error, id, curve = c('shcal20', 'intcal20', 'marine20',
               'intcal13', 'shcal13', 'marine13', 'intcal09'), eps = 1e-5) {
  curve <- read_curve(paste0(curve, '.14c'))
  cal_one <- function (age, error) {
    probs <- with(curve, {
      v = (age - rc_age)^2 / (error^2 + rc_err^2)
      exp(-(v / 2)) / sqrt(error^2 + rc_err^2)
    })
    probs <- normalize(probs)
    who <- probs >= eps
    pmf(value = curve[['year']][who], prob = probs[who])
  }
  if (length(age) != length(error)) stop('age and error lengths differ')
  n <- length(age)
  cal <- list()
  for (i in 1:n) cal[[i]] <- cal_one(age[i], error[i])
  if(!missing(id)) {
    if (length(age) != length(id)) stop('age and id lengths differ')
    names(cal) <- id
  }
  return(cal)
}

read_curve <- function(file) {
  path <- system.file('extdata', file, package = 'schron', mustWork = TRUE)
  curve <- read.table(path, sep = ',', col.names = c('cal', 'rc_age', 'rc_err',
    'delta', 'sigma'))
  year <- if(file %in% c('intcal20.14c', 'shcal20.14c', 'marine20.14c')) {
    seq(from = 55000, to = 0, by = -1)
  } else {
    seq(from = 50000, to = 0, by = -1)
  }
  toyear <- function (x) with(curve, approx(cal, x, xout = year)[['y']])
  rc_age <- toyear(curve[['rc_age']])
  rc_err <- toyear(curve[['rc_err']])
  list(year = year, rc_age = rc_age, rc_err = rc_err)
}