Nothing
rej_sampling = function(x,y, n = 1){
#' @noRd
#' @keywords internal
#'
#' @title rejection sampling
#'
#' @param x numeirc vector, strictly increasing.
#' @param y numeric vector, positive
#' @param n interger, number of samples generated
#'
#' @description
#' draws samples from a pdf given by approxfun(x,y)
#'
#' @returns numeric vector of length n
f = stats::approxfun(x,y)
f_max = max(y)
out = c()
while (length(out) < n){
x_draw = stats::runif(1, min = min(x), max = max(x))
y_draw = stats::runif(1, min = 0, max = f_max)
if (y_draw < f(x_draw)){
out = c(out, x_draw)
}
}
return(out)
}
crppp = function(x_min, x_max, rate = 1){
#' @noRd
#' @keywords internal
#'
#' @title Poisson point process
#'
#' @description
#' generates points of Poisson point process in the interval x_min, x_max
#'
#' @param x_min lower limit
#' @param x_max upper limit
#' @param rate numeric, rate of PPP
#'
#' @returns numeric vector of variable length
#'
n = stats::rpois(1, (x_max - x_min)* rate)
points = stats::runif(n, min = x_min, max = x_max)
return(points)
}
sed_rate_from_matrix = function(height, sedrate, matrix, mode = "deterministic", rate = 1, expand_domain = TRUE, transform = identity){
#' @export
#' @title make sed rate gen from matrix
#'
#' @param height vector of heights
#' @param sedrate vector of sed. rates x values
#' @param matrix matrix of sed rate y values. Must have as many columns as length(height) and as many rows as length(sedrate).
#' @param mode character, "deterministic" or "poisson". Determines at which stratigraphic heights the sed rate is determined. If "deterministic" this will be the heights in `height`, if "poisson" the heights where the sed rate is determined follows a poisson point process with rate specified by `rate`
#' @param rate numeric, rate of the Poisson point process determining frequency of sedimentation rate changes.
#' @param expand_domain should sedimentation rates be defined below/above the highest/lowest height in the section? If TRUE, the sed rate values are the values at the closest interpolated point, if FALSE it will be NA
#' @param transform a function, the identity function by default. How should the values of the (pseudo)pdf defined by the entries of `matrix` be transformed? Using this function allows to (nonlinearly) rescale the values in matrix to put more emphasis on higher/lower values
#'
#' @returns a function factory for usage with `sedrate_to_multiadm`
#'
#' @seealso [sedrate_to_multiadm()] for estimating sedimentation rates based on the outputs, [get_data_from_eTimeOpt()] for extracting data from the `eTimeOpt` function of the astrochron package.
#'
#' @description
#' Construct a sedimentation rate generator (function factory) from a matrix, e.g. one returned from `get_data_from_eTimeOpt`. This generator can be passed on to `sedrate_to_multiadm` to estimate age-depth models from it.
#' If mode is "deterministic", the generator evaluates the sedimentation rates at heights specified by `height`, if the mode is "poisson" it is evaluated at heights that are determined based on a poisson point process. At these heights, the value of the sedimentation rate is determined based on the (pseudo) pdf that is determined by the matrix values.
if (!all(dim(matrix) == c(length(sedrate), length(height)))){
stop("dimension mismatch. \"matrix\" must have length(height) columns and length(sedrate) rows")
}
rule = 1
if (expand_domain == TRUE){
rule = 2
}
if (mode == "poisson"){
f = function(){
x_max = max(height)
x_min = min(height)
interp_points = sort(c(x_min, x_max, crppp(x_min, x_max, rate)))
interp_heights = rep(NA, length(interp_points))
interp_vals = rep(NA, length(interp_points))
se = rep(NA, length(interp_points))
for (i in seq_along(interp_points)){
low_int = height[height <= interp_points[i]] #heights below poi
high_int = height[height >= interp_points[i]] # heights above poi
h_low = max(low_int) # highest point below poi
h_high = min(high_int) # lowest point above poi
low_ind = which(height ==h_low)
high_ind = which(height == h_high)
if (h_high == h_low){
lam = 0
} else {
lam = (interp_points[i] - h_low)/(h_high - h_low) #relative position of the poi between h_low and h_high
}
ppdf = lam * matrix[, low_ind] + (1-lam)* matrix[, high_ind]
sed_rate_vals = transform(ppdf)
sed_rate_val = rej_sampling(sedrate, sed_rate_vals)
se[i] = sed_rate_val
}
return(stats::approxfun(x = interp_points, y = se, ties = function(x) sample(x, 1), rule = rule)) # for ties, randomly select one sample
}
return(f)
}
if (mode == "deterministic"){
f = function(){
interp_points = sort(height)
se = rep(NA, length(interp_points))
for (i in seq_along(interp_points)){
sed_rate_vals = transform(matrix[, i])
se[i] = rej_sampling(sedrate, sed_rate_vals)
}
return(stats::approxfun(x = interp_points, y = se, rule = rule))
}
return(f)
}
stop("unknown mode, use either \"poisson\" or \"deterministic\"")
}
sed_rate_gen_from_bounds = function(h_l, s_l, h_u, s_u, rate = 1){
#' @export
#'
#' @title seg rate gen from upper/lower bounds
#'
#' @param h_l height values for lower bounds
#' @param s_l sed rate values for lower bounds
#' @param h_u height values for upper bounds
#' @param s_u sed rate values for upper bounds
#' @param rate rate of poisson point process
#'
#' @returns a function factory for usage with `sedrate_to_multiadm`
#'
#' @seealso [sedrate_to_multiadm()] for estimating age-depth models using the outputs, [sed_rate_from_matrix()] for other means of defining sedimentation rates
f = function(){
h_min = min(c(min(h_u), min(h_l)))
h_max = max(c(max(h_u), max(h_l)))
h = sort(unique(c(h_min, h_max, crppp(h_min, h_max, rate))))
s_min = stats::approx(h_l, s_l, xout = h, f = 2)$y
s_max = stats::approx(h_u, s_u, xout = h, f = 2)$y
sval = stats::runif(length(h), s_min, s_max)
return(stats::approxfun(h, sval, f=2))
}
return(f)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.