Nothing
#' Internal rust functions
#'
#' Internal rust functions
#' @details
#' These functions are not intended to be called by the user.
#' @name rust-internal
#' @keywords internal
NULL
# ==================================== box_cox ================================
#' @keywords internal
#' @rdname rust-internal
box_cox <- function (x, lambda = 1, gm = 1, lambda_tol = 1e-6) {
#
# Computes the Box-Cox transformation of a vector.
#
# Args:
# x : A numeric vector. (Positive) values to be Box-Cox
# transformed.
# lambda : A numeric scalar. Transformation parameter.
# gm : A numeric scalar. Optional scaling parameter.
# lambda_tol : A numeric scalar. For abs(lambda) < lambda_tol use
# a Taylor series expansion.
#
# Returns:
# A numeric vector. The transformed value
# (x^lambda - 1) / (lambda * gm ^ (lambda - 1))
#
if (abs(lambda) > lambda_tol) {
retval <- (x ^ lambda - 1) / lambda / gm ^ (lambda - 1)
} else {
i <- 0:3
retval <- sum(log(x) ^ (i + 1) * lambda ^ i / factorial(i + 1))
retval <- retval / gm ^ (lambda - 1)
}
retval
}
# ================================ box_cox_vec ================================
# Version of box_cox vectorized for lambda and gm.
#' @keywords internal
#' @rdname rust-internal
box_cox_vec <- Vectorize(box_cox, vectorize.args = c("lambda", "gm"))
# =============================== optim_box_cox ===============================
#' @keywords internal
#' @rdname rust-internal
optim_box_cox <- function(x, w, lambda_range = c(-3,3), start = NULL,
which_lam = 1:ncol(x)) {
#
# Finds the optimal value of the Box-Cox transformation parameter lambda,
# based values of a probability density function.
#
# Args:
# x : A numeric matrix. Values at which the density is evaluated.
# each row contains a combination of the ncol(x) variables.
# column numbers in which_lam must contain positive values.
# w : A numeric vector. Density values corresponding to each
# row of x (up to proportionality).
# lambda_range : A numeric vector (of length 2). Range of lambda values
# over which to search.
# start : A numeric vector. Optional starting value for lambda.
# which_lam : A numeric vector. Indicates which variables to Box-Cox
# transform.
#
# Returns: a list containing
# lambda : A numeric vector. The optimal value of lambda.
# gm : A numeric vector. Geometric mean of x weighted by w.
#
n_var <- ncol(x)
neg_loglik <- function(lambda_in, x, gm, w, which_lam) {
lambda <- rep(1, n_var)
lambda[which_lam] <- lambda_in
for (j in 1:n_var) {
x[, j] <- box_cox(x = x[, j], lambda = lambda[j], gm = gm[j])
}
(nrow(x) / 2) * log(det(stats::cov.wt(x, w, method = "ML")$cov))
}
w_n <- w / mean(w)
n_lam <- length(which_lam)
#
if (any(x[, which_lam] <= 0)) {
stop("All values must be > 0")
}
gm_w <- rep(1, n_var)
x_mat <- x[, which_lam, drop = FALSE]
gm_w[which_lam] <- apply(x_mat, 2, function(x) exp(mean(w_n * log(x))))
#
skew_wt <- function(x, w) {
vp <- function(p) mean(w^p)
v1 <- vp(1)
xbar <- mean(w * x) / v1
mp <- function(p) mean(w * (x - xbar) ^ p) / v1
mp(3) / mp(2)^(3 / 2)
}
if (is.null(start)) {
start <- 1 - apply(x, 2, skew_wt, w = w) / 2
}
if (n_lam == 1L) {
method <- "Brent"
} else {
method <- "L-BFGS-B"
}
lower <- lambda_range[1]
upper <- lambda_range[2]
start <- start[which_lam]
ret <- stats::optim(start, neg_loglik, method = method, x = x, gm = gm_w,
w = w_n, lower = lower, upper = upper,
which_lam = which_lam)
lambda <- rep(1L, n_var)
lambda[which_lam] <- ret$par
ret$par <- lambda
if (ret$convergence != 0) {
warning(paste("Convergence failure: return code =", ret$convergence))
}
list(lambda = ret$par, gm = gm_w)
}
# ================================ n_grid_fun =================================
#' @keywords internal
#' @rdname rust-internal
n_grid_fn <- function(d) {
# Sets a default values value of n_grid.
#
# Args:
# d : A numeric scalar. Dimension of the target density.
#
# Returns:
# A numeric scalar. The value of n_grid.
#
return(ceiling(2501 ^ (1 / d)))
}
# =============================== init_ps_calc ================================
#' @keywords internal
#' @rdname rust-internal
init_psi_calc <- function(phi_to_psi, phi, lambda, gm, w, which_lam){
# Estimates the mode and standard deviation of the Box-Cox transformed
# target density.
#
# Args:
# phi_to_psi : A function. The Box-Cox transformation function.
# phi : A numeric matrix. n_grid by d matrix of values of phi
# lambda : A numeric vector. The value of the Box-Cox transformation
# parameter lambda.
# gm : A numeric vector. The value of the Box-Cox scale
# parameter gm.
# w : A numeric vector. The values of the target density at each
# value of the variables in phi.
# which_lam : A numeric vector. Indicates which variables to Box-Cox
# transform.
#
# Returns: A list containing
# init_psi : estimate of the mode of the Box-Cox transformed target
# density.
# sd_psi : estimate of the standard deviation of the Box-Cox transformed
# target density.
#
psi <- phi_to_psi(phi)
d <- ncol(phi)
log_jac <- matrix(0, nrow(phi), d)
for (j in which_lam) {
log_jac[, j] <- (lambda[j] - 1) * log(phi[, j] / gm[j])
}
log_jac <- rowSums(log_jac)
w_psi <- w * exp(-log_jac)
temp_cov <- stats::cov.wt(psi, w_psi, method="ML")
init_psi <- as.numeric(psi[which.max(w_psi), ])
sd_psi <- sqrt(diag(temp_cov$cov))
list(init_psi = init_psi, sd_psi = sd_psi)
}
# ============================= log_gpd_mdi_prior =============================
#' @keywords internal
#' @rdname rust-internal
log_gpd_mdi_prior <- function(pars) {
# Generalized Pareto MDI prior log-density
#
# Calculates the log-density of a Maximal Data Information (MDI) prior
# for generalized Pareto parameters. The prior density is truncated,
# (to xi >= -1), to produce a posterior density that is proper.
#
# Args:
# pars : A numeric vector containing the values of the generalized Pareto
# parameters sigma and xi.
#
# Returns:
# A numeric scalar. The value of the prior log-density.
#
if (pars[1] <= 0 | pars[2] < -1) {
return(-Inf)
}
return(-log(pars[1]) - pars[2] - 1)
}
# ================================ gpd_loglik =================================
#' @keywords internal
#' @rdname rust-internal
gpd_loglik <- function(pars, gpd_data, m, xm, sum_gp) {
# Generalized Pareto log-likelihood
#
# Calculates the log-likelihood for a random sample from a Generalized Pareto.
#
# Args:
# pars : A numeric vector containing the values of the generalized
# Pareto parameters sigma and xi.
# gpd_data : A numeric vector containing positive sample values.
# m : A numeric scalar. The sample size: the length of gpd_data.
# xm : A numeric scalar. The sample maximum.
# sum_gp : The sum of the sample values.
#
# Returns:
# A numeric scalar. The value of the log-likelihood.
#
if (pars[1] <= 0 | pars[2] <= -pars[1] / xm)
return(-Inf)
sdat <- gpd_data / pars[1]
zz <- 1 + pars[2] * sdat
if (abs(pars[2]) > 1e-6) {
val <- -m * log(pars[1]) - (1 + 1 / pars[2]) * sum(log(zz))
} else {
i <- 1:4
g_fun <- function(x) {
t1 <- x ^ i
t2 <- (i * x - i - 1)
sum((-1) ^ i * t1 * t2 * pars[2] ^ i / i / (i + 1))
}
val <- -m * log(pars[1]) - sum_gp / pars[1] - sum(sapply(sdat, g_fun))
}
return(val)
}
# ================================== gpd_mle ==================================
#' @keywords internal
#' @rdname rust-internal
gpd_mle <- function(gpd_data) {
# Maximum likelihood estimation for the generalized Pareto distribution
#
# Performs maximum likelihood estimation for the generalized Pareto
# distribution. Uses the function \code{gpdmle} associated with
# Grimshaw (1993), which returns MLEs of sigma and k = - \xi.
#
# Grimshaw, S. D. (1993) Computing Maximum Likelihood Estimates
# for the Generalized Pareto Distribution. Technometrics, 35(2), 185-191.
# and Computing (1991) 1, 129-133. https://doi.org/10.1007/BF01889987.
#
# Args:
# gpd_data : A numeric vector containing positive values, assumed to be a
# random sample from a generalized Pareto distribution.
#
# Returns:
# A list with components
# mle : A numeric vector. MLEs of GP parameters sigma and xi.
# nllh : A numeric scalar. The negated log-likelihood at the MLE.
#
temp <- list()
# Use revdbayes function if revdbayes is available, otherwise use fallback
if (requireNamespace("revdbayes", quietly = TRUE)) {
# Call Grimshaw (1993) function, note: k is -xi, a is sigma
pjn <- revdbayes::grimshaw_gp_mle(gpd_data)
temp$mle <- c(pjn$a, -pjn$k) # mle for (sigma,xi)
} else {
temp <- fallback_gp_mle(init = c(mean(gpd_data), 0), gpd_data = gpd_data,
m = length(gpd_data), xm = max(gpd_data),
sum_gp = sum(gpd_data))
}
sc <- rep(temp$mle[1], length(gpd_data))
xi <- temp$mle[2]
temp$nllh <- sum(log(sc)) + sum(log1p(xi * gpd_data / sc) * (1 / xi + 1))
return(temp)
}
# ============================== fallback_gp_mle ==============================
#' @keywords internal
#' @rdname rust-internal
fallback_gp_mle <- function(init, ...){
x <- stats::optim(init, gpd_loglik, ..., control = list(fnscale = -1),
hessian = FALSE)
temp <- list()
temp$mle <- x$par
temp$nllh <- -x$value
return(temp)
}
# =========================== gpd_obs_info ===========================
#' @keywords internal
#' @rdname rust-internal
gpd_obs_info <- function(gpd_pars, y, eps = 1e-5, m = 3) {
# Observed information for the generalized Pareto distribution
#
# Calculates the observed information matrix for a random sample \code{y}
# from the generalized Pareto distribution, i.e. the negated Hessian matrix
# of the generalized Pareto log-likelihood, evaluated at \code{gp_pars}.
#
# Args:
# gp_pars : A numeric vector. Parameters sigma and xi of the
# generalized Pareto distribution.
# y : A numeric vector. A sample of positive data values.
# eps : A (small, positive) numeric scalar. If abs(xi) is smaller than
# eps then we approximate the [2, 2] element of the information
# matrix using a Taylor series approximation.
# m : A non-negative integer. The order of the polynomial used to
# make this approximation.
# Returns:
# A 2 by 2 numeric matrix. The observed information matrix.
#
if (eps <= 0) {
stop("'eps' must be positive")
}
if (m < 0) {
stop("'m' must be non-negative")
}
# sigma
s <- gpd_pars[1]
# xi
x <- gpd_pars[2]
i <- matrix(NA, 2, 2)
i[1, 1] <- -sum((1 - (1 + x) * y * (2 * s + x * y) / (s + x * y) ^ 2) / s ^ 2)
i[1, 2] <- i[2, 1] <- -sum(y * (1 - y / s) / (1 + x * y / s) ^ 2 / s ^ 2)
# Direct calculation of i22 is unreliable for x close to zero.
# If abs(x) < eps then we expand the problematic terms (all but t4 below)
# in powers of xi up to xi ^ m. The terms in 1/z and 1/z^2 cancel.
z <- x / s
t0 <- 1 + z * y
t4 <- y ^ 2 / t0 ^ 2
if (any(t0 <= 0)) {
stop("The log-likelihood is 0 for this combination of data and parameters")
}
if (abs(x) < eps) {
j <- 0:m
zy <- z * y
sum_fn <- function(zy) {
return(sum((-1) ^ j * (j ^ 2 + 3 * j + 2) * zy ^ j / (j + 3)))
}
tsum <- vapply(zy, sum_fn, 0.0)
i[2, 2] <- sum(y ^ 3 * tsum / s ^ 3 - t4 / s ^ 2)
} else {
t1 <- 2 * log(t0) / z ^ 3
t2 <- 2 * y / (z ^ 2 * t0)
t3 <- y ^ 2 / (z * t0 ^ 2)
i[2, 2] <- sum((t1 - t2 - t3) / s ^ 3 - t4 / s ^ 2)
}
dimnames(i) <- list(c("sigma[u]", "xi"), c("sigma[u]", "xi"))
return(i)
}
# =================================== find_a ==================================
#' @keywords internal
#' @rdname rust-internal
find_a <- function(neg_logf_rho, init_psi, d, r, lower, upper, algor,
method, control, shoof, ...) {
#
# Finds the value of a(r).
#
# Args:
# neg_logf_rho : A function. Negated target log-density function.
# init_psi : A numeric scalar. Initial value of psi.
# d : A numeric scalar. Dimension of f.
# r : A numeric scalar. Parameter of generalized
# ratio-of-uniforms.
# lower : A numeric vector. Lower bounds on the arguments of logf.
# upper : A numeric vector. Upper bounds on the arguments of logf.
# algor : A character scalar. Algorithm ("optim" or "nlminb").
# method : A character scalar. Only relevant if algorithm = "optim".
# control : A numeric list. Control arguments to algor.
#
# Returns: a list containing
# the standard returns from optim or nlminb
# hessian: the estimated hessian of neg_logf_rho/(d*r+1) at its minimum.
#
big_val <- Inf
big_finite_val <- 10 ^ 10
#
# Function to minimize to find a(r).
# Use a weird argument name (._psi) to avoid a potential argument-matching
# problem when nlminb() is used.
a_obj <- function(._psi, ...) {
# Avoid possible problems with nlminb calling function with NaNs.
# Indicative of solution on boundary, which is OK in the current context.
# See https://stat_ethz.ch/pipermail/r-help/2015-May/428488.html
if (any(is.na(._psi))) return(big_val)
neg_logf_rho(._psi, ...) / (d * r + 1)
}
# Function for L-BFGS-B and Brent, which don't like Inf or NA
a_obj_no_inf <- function(._psi, ...) {
check <- neg_logf_rho(._psi, ...) / (d * r + 1)
if (is.infinite(check)) {
check <- big_finite_val
}
check
}
#
if (algor == "optim") {
if (method %in% c("L-BFGS-B","Brent")) {
temp <- stats::optim(par = init_psi, fn = a_obj_no_inf, ...,
control = control, hessian = FALSE, method = method,
lower = lower, upper = upper)
} else {
temp <- stats::optim(par = init_psi, fn = a_obj, ..., control = control,
hessian = FALSE, method = method)
# Sometimes Nelder-Mead fails if the initial estimate is too good.
# ... so avoid non-zero convergence indicator by using L-BFGS-B instead.
if (temp$convergence > 0) {
# Start a little away from the optimum, to avoid erroneous
# convergence warnings, using init_psi as a benchmark
# If init_psi = temp$par then multiply temp$par by 1 - shoof
if (sum(abs(init_psi - temp$par)) > .Machine$double.eps) {
new_start <- shoof * init_psi + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
temp <- stats::optim(par = new_start, fn = a_obj_no_inf, ...,
control = control, hessian = FALSE,
method = "L-BFGS-B", lower = lower, upper = upper)
}
# In some cases optim with method = "L-BFGS-B" may reach its iteration
# limit without the convergence criteria being satisfied. Then try
# nlminb as a further check, but don't use the control argument in
# case of conflict between optim() and nlminb().
if (temp$convergence > 0) {
temp <- stats::nlminb(start = new_start, objective = a_obj, ...,
lower = lower, upper = upper)
}
}
# Try to calculate Hessian, but avoid problems if an error is produced.
# An error may occur if the MAP estimate is very close to a parameter
# boundary.
temp$hessian <- try(stats::optimHess(par = temp$par, fn = a_obj, ...),
silent = TRUE)
return(temp)
}
# If we get to here we are using nlminb() ...
temp <- stats::nlminb(start = init_psi, objective = a_obj, ...,
lower = lower, upper = upper, control = control)
# Sometimes nlminb isn't sure that it has found the minimum when in fact
# it has. Try to check this, and avoid a non-zero convergence indicator
# by using optim with method="L-BFGS-B", again starting from new_start,
# but don't use the control argument in case of conflict between
# optim() and nlminb().
if (temp$convergence > 0) {
if (sum(abs(init_psi - temp$par)) > .Machine$double.eps) {
new_start <- shoof * init_psi + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
temp <- stats::optim(par = new_start, fn = a_obj_no_inf, ...,
hessian = FALSE, method = "L-BFGS-B", lower = lower,
upper = upper)
}
# Try to calculate Hessian, but avoid problems if an error is produced.
# An error may occur if the MAP estimate is very close to a parameter
# boundary.
temp$hessian <- try(stats::optimHess(par = temp$par, fn = a_obj, ...),
silent = TRUE)
return(temp)
}
# ================================== find_bs ==================================
#' @keywords internal
#' @rdname rust-internal
find_bs <- function(f_rho, d, r, lower, upper, f_mode, ep, vals, conv, algor,
method, control, shoof, ...) {
# Finds the values of b-(r) and b+(r).
#
# Args:
# f_rho : A function. Target probability density function.
# d : A numeric scalar. Dimension of f.
# r : A numeric scalar. Parameter of generalized
# ratio-of-uniforms.
# lower : A numeric vector. Lower bounds on the arguments of logf.
# upper : A numeric vector. Upper bounds on the arguments of logf.
# f_mode : A numeric scalar. The estimated mode of the target
# log-density logf.
# ep : A numeric scalar. Controls initial estimates for
# optimisations to find the b-bounding box parameters.
# The default (ep=0) corresponds to starting at the mode of
# logf small positive values of ep move the constrained
# variable slightly away from the mode in the correct
# direction. If ep is negative its absolute value is used,
# with no warning given.
# vals : A numeric matrix. Will contain the values of the
# variables at which the ru box dimensions occur.
# Row 1 already contains the values for a(r).
# conv : A numeric scalar. Will contain the convergence
# indicators returned by the optimisation algorithms.
# Row 1 already contains the values for a(r).
# algor : A character scalar. Algorithm ("optim" or "nlminb").
# method : A character scalar. Only relevant if algorithm = "optim".
# control : A numeric list. Control arguments to algor.
#
# Returns: a list containing
# l_box : A numeric vector. Values of biminus(r), i = 1, ...d.
# u_box : A numeric vector. Values of biplus(r), i = 1, ...d.
# vals : as described above in Args.
# conv : as described above in Args.
#
big_val <- Inf
big_finite_val <- 10^10
#
# Functions to minimize to find biminus(r) and biplus(s), i = 1, ...,d.
# Use a weird argument name (._rho) to avoid a potential argument-matching
# problem when nlminb() is used.
lower_box <- function(._rho, j, ...) {
# Avoid possible problems with nlminb calling function with NaNs.
# Indicative of solution on boundary, which is OK in the current context.
# See https://stat_ethz.ch/pipermail/r-help/2015-May/428488.html
if (any(is.na(._rho))) return(big_val)
if (._rho[j] == 0L) return(0L)
if (._rho[j] > 0L) return(big_val)
if (f_rho(._rho, ...) == 0L) return(big_val)
._rho[j] * f_rho(._rho, ...) ^ (r / (d * r + 1))
}
upper_box <- function(._rho, j, ...) {
if (any(is.na(._rho))) return(big_val)
if (._rho[j] == 0) return(0)
if (._rho[j] < 0) return(big_val)
if (f_rho(._rho, ...) == 0) return(big_val)
-._rho[j] * f_rho(._rho, ...) ^ (r / (d * r + 1))
}
# Functions for L-BFGS-B and Brent, which don't like Inf or NA
lower_box_no_inf <- function(._rho, j, ...) {
if (any(is.na(._rho))) return(big_finite_val)
if (._rho[j] == 0) return(0)
if (._rho[j] > 0) return(big_finite_val)
if (f_rho(._rho, ...) == 0) return(big_finite_val)
check <- ._rho[j] * f_rho(._rho, ...) ^ (r / (d * r + 1))
if (is.infinite(check)) check <- big_finite_val
check
}
upper_box_no_inf <- function(._rho, j, ...) {
if (any(is.na(._rho))) return(big_finite_val)
if (._rho[j] == 0) return(0)
if (._rho[j] < 0) return(big_finite_val)
if (f_rho(._rho, ...) == 0) return(big_finite_val)
check <- -._rho[j] * f_rho(._rho, ...) ^ (r / (d * r + 1))
if (is.infinite(check)) check <- big_finite_val
check
}
l_box <- u_box <- NULL
zeros <- rep(0,d)
#
# Find biminus(r) and biplus(s), i = 1, ...,d.
#
for (j in 1:d) {
#
# Find biminus(r) ----------
#
rho_init <- zeros
rho_init[j] <- -ep
t_upper <- upper - f_mode
t_upper[j] <- 0
if (algor == "nlminb") {
temp <- stats::nlminb(start = rho_init, objective = lower_box, j = j,
..., upper = t_upper, lower = lower - f_mode,
control = control)
l_box[j] <- temp$objective
# Sometimes nlminb isn't sure that it has found the minimum when in fact
# it has. Try to check this, and avoid a non-zero convergence indicator
# by using optim with method="L-BFGS-B", starting from nlminb's solution.
if (temp$convergence > 0) {
if (sum(abs(rho_init - temp$par)) > .Machine$double.eps) {
new_start <- shoof * rho_init + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
temp <- stats::optim(par = new_start, fn = lower_box_no_inf, j = j, ...,
hessian = FALSE, method = "L-BFGS-B",
upper = t_upper, lower = lower - f_mode)
l_box[j] <- temp$value
}
}
if (algor == "optim") {
if (method == "L-BFGS-B" | method == "Brent") {
temp <- stats::optim(par = rho_init, fn = lower_box_no_inf, j = j, ...,
upper = t_upper, lower = lower - f_mode,
control = control, method = method,
hessian = FALSE)
l_box[j] <- temp$value
} else {
temp <- stats::optim(par = rho_init, fn = lower_box, j = j, ...,
control = control, method = method,
hessian = FALSE)
l_box[j] <- temp$value
# Sometimes Nelder-Mead fails if the initial estimate is too good.
# ... so avoid non-zero convergence indicator using L-BFGS-B instead.
if (temp$convergence > 0) {
if (sum(abs(rho_init - temp$par)) > .Machine$double.eps) {
new_start <- shoof * rho_init + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
temp <- stats::optim(par = new_start, fn = lower_box_no_inf, j = j,
..., control = control, method = "L-BFGS-B",
hessian = FALSE,
upper = t_upper, lower = lower - f_mode)
l_box[j] <- temp$value
}
# Check using nlminb() if optim's iteration limit is reached.
if (temp$convergence == 1) {
if (sum(abs(rho_init - temp$par)) > .Machine$double.eps) {
new_start <- shoof * rho_init + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
temp <- stats::nlminb(start = new_start, objective = lower_box,
j = j, ..., upper = t_upper,
lower = lower - f_mode)
l_box[j] <- temp$objective
}
}
}
vals[j+1, ] <- temp$par
conv[j+1] <- temp$convergence
#
# Find biplus(r) --------------
#
rho_init <- zeros
rho_init[j] <- ep
t_lower <- lower - f_mode
t_lower[j] <- 0
if (algor == "nlminb") {
temp <- stats::nlminb(start = rho_init, objective = upper_box, j = j,
..., lower = t_lower, upper = upper - f_mode,
control = control)
u_box[j] <- -temp$objective
# Sometimes nlminb isn't sure that it has found the minimum when in fact
# it has. Try to check this, and avoid a non-zero convergence indicator
# by using optim with method="L-BFGS-B", starting from nlminb's solution.
if (temp$convergence > 0) {
if (sum(abs(rho_init - temp$par)) > .Machine$double.eps) {
new_start <- shoof * rho_init + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
temp <- stats::optim(par = new_start, fn = upper_box_no_inf, j = j,
..., hessian = FALSE, method = "L-BFGS-B",
lower = t_lower, upper = upper - f_mode)
u_box[j] <- -temp$value
}
}
if (algor == "optim") {
if (method == "L-BFGS-B" | method == "Brent") {
temp <- stats::optim(par = rho_init, fn = upper_box_no_inf, j = j, ...,
lower = t_lower, upper = upper - f_mode,
control = control, method = method)
u_box[j] <- -temp$value
} else {
temp <- stats::optim(par = rho_init, fn = upper_box, j = j, ...,
control = control, method = method)
u_box[j] <- -temp$value
# Sometimes Nelder-Mead fails if the initial estimate is too good.
# ... so avoid non-zero convergence indicator using L-BFGS-B instead.
if (temp$convergence > 0) {
if (sum(abs(rho_init - temp$par)) > .Machine$double.eps) {
new_start <- shoof * rho_init + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
temp <- stats::optim(par = new_start, fn = upper_box_no_inf, j = j,
..., control = control, method = "L-BFGS-B",
lower = t_lower, upper = upper - f_mode)
u_box[j] <- -temp$value
}
# Check using nlminb() if optim's iteration limit is reached.
if (temp$convergence == 1) {
if (sum(abs(rho_init - temp$par)) > .Machine$double.eps) {
new_start <- shoof * rho_init + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
temp <- stats::nlminb(start = new_start, objective = upper_box,
j = j, ..., lower = t_lower,
upper = upper - f_mode)
u_box[j] <- -temp$objective
}
}
}
vals[j+d+1, ] <- temp$par
conv[j+d+1] <- temp$convergence
}
return(list(l_box = l_box, u_box = u_box, vals = vals, conv = conv))
}
# ================================ cpp_find_a =================================
#' @keywords internal
#' @rdname rust-internal
cpp_find_a <- function(init_psi, lower, upper, algor, method, control,
a_obj_fun, ru_args, shoof) {
#
# Finds the value of a(r).
#
# Args:
# init_psi : A numeric scalar. Initial value of psi.
# lower : A numeric vector. Lower bounds on the arguments of logf.
# upper : A numeric vector. Upper bounds on the arguments of logf.
# algor : A character scalar. Algorithm ("optim" or "nlminb").
# method : A character scalar. Only relevant if algorithm = "optim".
# control : A numeric list. Control arguments to algor.
# a_obj_fun : The function to be minimized to find a(r).
# ru_args : A numeric list, containing:
# d : A numeric scalar. Dimension of f.
# r : A numeric scalar. Generalized ratio-of-uniforms parameter.
# psi_mode : A numeric vector. Latterly this will contain the estimated
# mode of the target density after transformation but prior
# to mode relation. Equal to rep(0, d) at this stage.
# rot_mat : A numeric matrix. Latterly this will contain the rotation
# matrix. Equal to diag(d) at this stage.
# hscale : A numeric scalar. Scales the target log-density.
# Equal to logf evaluated at init at this stage.
# which_lam : A vector of integers indication which components of lambda
# are Box-Cox transformed. Only present if trans = "BC".
# lambda : Box-Cox transformation parameters.
# Only present if trans = "BC".
# gm : Box-Cox scale parameters. Only present if trans = "BC".
# con : lambda * gm ^(lambda - 1). Only present if trans = "BC".
# logf : A pointer to the (original) target log-density function.
# pars : A numeric list of additional arguments to logf.
#
# Returns: a list containing
# the standard returns from optim or nlminb
# hessian: the estimated hessian of -cpp_logf_rho/(d*r+1) at its minimum.
#
big_val <- 10 ^ 10
#
if (algor == "optim") {
if (method == "L-BFGS-B" | method == "Brent") {
add_args <- list(par = init_psi, fn = a_obj_fun, method = method,
control = control, lower = lower, upper = upper,
big_val = big_val)
temp <- do.call(stats::optim, c(ru_args, add_args))
} else {
add_args <- list(par = init_psi, fn = a_obj_fun, method = method,
control = control, big_val = Inf)
temp <- do.call(stats::optim, c(ru_args, add_args))
# Sometimes Nelder-Mead fails if the initial estimate is too good.
# ... so avoid non-zero convergence indicator using L-BFGS-B instead.
if (temp$convergence > 0) {
# Start a little away from the optimum, to avoid erroneous
# convergence warnings, using init_psi as a benchmark
# If init_psi = temp$par then multiply temp$par by 1 - shoof
if (sum(abs(init_psi - temp$par)) > .Machine$double.eps) {
new_start <- shoof * init_psi + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
add_args <- list(par = new_start, fn = a_obj_fun, method = "L-BFGS-B",
control = control, big_val = big_val,
lower = lower, upper = upper)
temp <- do.call(stats::optim, c(ru_args, add_args))
}
# In some cases optim with method = "L-BFGS-B" may reach its iteration
# limit without the convergence criteria being satisfied. Then try
# nlminb as a further check, but don't use the control argument in
# case of conflict between optim() and nlminb().
if (temp$convergence > 0) {
add_args <- list(start = new_start, objective = a_obj_fun,
lower = lower, upper = upper, big_val = Inf)
temp <- do.call(stats::nlminb, c(ru_args, add_args))
}
}
} else {
add_args <- list(start = init_psi, objective = a_obj_fun,
control = control, lower = lower, upper = upper,
big_val = Inf)
temp <- do.call(stats::nlminb, c(ru_args, add_args))
# Sometimes nlminb isn't sure that it has found the minimum when in fact
# it has. Try to check this, and avoid a non-zero convergence indicator
# by using optim with method="L-BFGS-B", again starting from new_start.
if (temp$convergence > 0) {
if (sum(abs(init_psi - temp$par)) > .Machine$double.eps) {
new_start <- shoof * init_psi + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
add_args <- list(par = new_start, fn = a_obj_fun, hessian = FALSE,
method = "L-BFGS-B", big_val = big_val,
lower = lower, upper = upper)
temp <- do.call(stats::optim, c(ru_args, add_args))
}
}
# Try to calculate Hessian, but avoid problems if an error is produced.
# An error may occur if the MAP estimate is very close to a parameter
# boundary.
add_args <- list(par = temp$par, fn = a_obj_fun, big_val = Inf)
temp$hessian <- try(do.call(stats::optimHess, c(ru_args, add_args)),
silent = TRUE)
return(temp)
}
# ================================ cpp_find_bs ================================
#' @keywords internal
#' @rdname rust-internal
cpp_find_bs <- function(lower, upper, ep, vals, conv, algor, method,
control, lower_box_fun, upper_box_fun, ru_args,
shoof) {
# Finds the values of b-(r) and b+(r).
#
# Args:
# lower : A numeric vector. Lower bounds on the arguments of logf.
# upper : A numeric vector. Upper bounds on the arguments of logf.
# ep : A numeric scalar. Controls initial estimates for
# optimisations to find the b-bounding box parameters.
# The default (ep=0) corresponds to starting at the mode of
# logf small positive values of ep move the constrained
# variable slightly away from the mode in the correct
# direction. If ep is negative its absolute value is used,
# with no warning given.
# vals : A numeric matrix. Will contain the values of the
# variables at which the ru box dimensions occur.
# Row 1 already contains the values for a(r).
# conv : A numeric scalar. Will contain the convergence
# indicators returned by the optimisation algorithms.
# Row 1 already contains the values for a(r).
# algor : A character scalar. Algorithm ("optim" or "nlminb").
# method : A character scalar. Only relevant if algor = "optim".
# control : A numeric list. Control arguments to algor.
# lower_box_fun : The function to be minimized to find b-(r).
# upper_box_fun : The function to be minimized to find b+(r).
# ru_args : A numeric list, containing:
# d : A numeric scalar. Dimension of f.
# r : A numeric scalar. Generalized ratio-of-uniforms parameter.
# psi_mode : A numeric vector. The estimated mode of the target
# density after transformation but prior to mode relocation.
# rot_mat : A numeric matrix. Rotation matrix (equal to the identity
# matrix if rotate = FALSE).
# hscale : A numeric scalar. Scales the target log-density.
# which_lam : A vector of integers indication which components of lambda
# are Box-Cox transformed. Only present if trans = "BC".
# lambda : Box-Cox transformation parameters.
# Only present if trans = "BC".
# gm : Box-Cox scale parameters. Only present if trans = "BC".
# con : lambda * gm ^(lambda - 1). Only present if trans = "BC".
# logf : A pointer to the (original) target log-density function.
# pars : A numeric list of additional arguments to logf.
#
# Returns: a list containing
# l_box : A numeric vector. Values of biminus(r), i = 1, ...d.
# u_box : A numeric vector. Values of biplus(r), i = 1, ...d.
# vals : as described above in Args.
# conv : as described above in Args.
#
big_val <- 10 ^ 10
f_mode <- ru_args$psi_mode
d <- ru_args$d
#
l_box <- u_box <- NULL
zeros <- rep(0, d)
#
# Find biminus(r) and biplus(s), i = 1, ...,d.
#
for (j in 1:d) {
#
# Find biminus(r) ----------
#
rho_init <- zeros
rho_init[j] <- -ep
t_upper <- upper - f_mode
t_upper[j] <- 0
if (algor == "nlminb") {
add_args <- list(start = rho_init, objective = lower_box_fun,
upper = t_upper, lower = lower - f_mode, j = j - 1,
control = control, big_val = Inf)
temp <- do.call(stats::nlminb, c(ru_args, add_args))
l_box[j] <- temp$objective
# Sometimes nlminb isn't sure that it has found the minimum when in fact
# it has. Try to check this, and avoid a non-zero convergence indicator
# by using optim with method="L-BFGS-B", starting from nlminb's solution.
if (temp$convergence > 0) {
if (sum(abs(rho_init - temp$par)) > .Machine$double.eps) {
new_start <- shoof * rho_init + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
add_args <- list(par = new_start, fn = lower_box_fun, j = j - 1,
method = "L-BFGS-B", big_val = big_val,
upper = t_upper, lower = lower - f_mode)
temp <- do.call(stats::optim, c(ru_args, add_args))
l_box[j] <- temp$value
}
}
if (algor == "optim") {
# L-BFGS-B and Brent don't like Inf or NA
if (method == "L-BFGS-B" | method == "Brent") {
add_args <- list(par = rho_init, fn = lower_box_fun, upper = t_upper,
lower = lower - f_mode, j = j - 1, control = control,
method = method, big_val = big_val)
temp <- do.call(stats::optim, c(ru_args, add_args))
l_box[j] <- temp$value
} else {
add_args <- list(par = rho_init, fn = lower_box_fun, j = j - 1,
control = control, method = method, big_val = Inf)
temp <- do.call(stats::optim, c(ru_args, add_args))
l_box[j] <- temp$value
# Sometimes Nelder-Mead fails if the initial estimate is too good.
# ... so avoid non-zero convergence indicator using L-BFGS-B instead.
if (temp$convergence > 0) {
if (sum(abs(rho_init - temp$par)) > .Machine$double.eps) {
new_start <- shoof * rho_init + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
add_args <- list(par = new_start, fn = lower_box_fun, j = j - 1,
control = control, method = "L-BFGS-B",
big_val = big_val, upper = t_upper,
lower = lower - f_mode)
temp <- do.call(stats::optim, c(ru_args, add_args))
l_box[j] <- temp$value
}
# Check using nlminb() if optim's iteration limit is reached.
if (temp$convergence == 1) {
if (sum(abs(rho_init - temp$par)) > .Machine$double.eps) {
new_start <- shoof * rho_init + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
add_args <- list(start = new_start, objective = lower_box_fun,
upper = t_upper, lower = lower - f_mode, j = j - 1,
big_val = Inf)
temp <- do.call(stats::nlminb, c(ru_args, add_args))
l_box[j] <- temp$objective
}
}
}
vals[j+1, ] <- temp$par
conv[j+1] <- temp$convergence
#
# Find biplus(r) --------------
#
rho_init <- zeros
rho_init[j] <- ep
t_lower <- lower - f_mode
t_lower[j] <- 0
if (algor == "nlminb") {
add_args <- list(start = rho_init, objective = upper_box_fun,
lower = t_lower, upper = upper - f_mode, j = j - 1,
control = control, big_val = Inf)
temp <- do.call(stats::nlminb, c(ru_args, add_args))
u_box[j] <- -temp$objective
# Sometimes nlminb isn't sure that it has found the minimum when in fact
# it has. Try to check this, and avoid a non-zero convergence indicator
# by using optim with method="L-BFGS-B", starting from nlminb's solution.
if (temp$convergence > 0) {
if (sum(abs(rho_init - temp$par)) > .Machine$double.eps) {
new_start <- shoof * rho_init + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
add_args <- list(par = new_start, fn = upper_box_fun, j = j - 1,
method = "L-BFGS-B", big_val = big_val,
lower = t_lower, upper = upper - f_mode)
temp <- do.call(stats::optim, c(ru_args, add_args))
u_box[j] <- -temp$value
}
}
if (algor == "optim") {
# L-BFGS-B and Brent don't like Inf or NA
if (method == "L-BFGS-B" | method == "Brent") {
add_args <- list(par = rho_init, fn = upper_box_fun,
lower = t_lower, upper = upper - f_mode, j = j - 1,
control = control, method = method, big_val = big_val)
temp <- do.call(stats::optim, c(ru_args, add_args))
u_box[j] <- -temp$value
} else {
add_args <- list(par = rho_init, fn = upper_box_fun, j = j - 1,
control = control, method = method, big_val = Inf)
temp <- do.call(stats::optim, c(ru_args, add_args))
u_box[j] <- -temp$value
# Sometimes Nelder-Mead fails if the initial estimate is too good.
# ... so avoid non-zero convergence indicator using L-BFGS-B instead.
if (temp$convergence > 0) {
if (sum(abs(rho_init - temp$par)) > .Machine$double.eps) {
new_start <- shoof * rho_init + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
add_args <- list(par = new_start, fn = upper_box_fun, j = j - 1,
control = control, method = "L-BFGS-B",
big_val = big_val, lower = t_lower,
upper = upper - f_mode)
temp <- do.call(stats::optim, c(ru_args, add_args))
u_box[j] <- -temp$value
}
# Check using nlminb() if optim's iteration limit is reached.
if (temp$convergence == 1) {
if (sum(abs(rho_init - temp$par)) > .Machine$double.eps) {
new_start <- shoof * rho_init + (1 - shoof) * temp$par
} else {
new_start <- temp$par * (1 - shoof)
}
add_args <- list(start = new_start, objective = upper_box_fun,
lower = t_lower, upper = upper - f_mode, j = j - 1,
big_val = Inf)
temp <- do.call(stats::nlminb, c(ru_args, add_args))
u_box[j] <- -temp$objective
}
}
}
vals[j+d+1, ] <- temp$par
conv[j+d+1] <- temp$convergence
}
return(list(l_box = l_box, u_box = u_box, vals = vals, conv = conv))
}
# ================================= wecdf =====================================
wecdf <- function (x, weights = rep(1, length(x))) {
# Sort x and reorder the weights
w <- weights[order(x)]
x <- sort(x)
# The sum of the weights for each unique value in xs
ws <- tapply(w, x, sum)
vals <- unique(x)
rval <- stats::approxfun(vals, cumsum(ws) / sum(ws), method = "constant",
yleft = 0, yright = 1, f = 0, ties = "ordered")
class(rval) <- c("ecdf", "stepfun", class(rval))
assign("nobs", length(x), envir = environment(rval))
attr(rval, "call") <- sys.call()
rval
}
# ============================ 3^d factorial design ======================== #
#' @keywords internal
#' @rdname rust-internal
fac3 <- function(d) {
# Creates a design matrix for a 3^d full factorial design
#
# Args:
# d : a positive integer scalar
#
# Returns:
# a 3^d by d integer matrix containing -1, 0, 1 for the levels of the d
# factors. The first column varies fastest, the second column second
# fastest etc.
rep_fn <- function(x) {
rep(-1L:1L, times = 3L ^ x, each = 3L ^ (d - 1L - x))
}
return(vapply((d - 1L):0, rep_fn, numeric(3L ^ d)))
}
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.