R/box_cox_functions_rcpp.R

In rust: Ratio-of-Uniforms Simulation with Transformation

# ========================= find_lambda_one_d_rcpp =========================

#' Selecting the Box-Cox parameter in the 1D case using Rcpp
#'
#' Finds a value of the Box-Cox transformation parameter lambda for which
#' the (positive univariate) random variable with log-density
#' \eqn{\log f}{log f} has a density closer to that of a Gaussian random
#' variable. Works by estimating a set of quantiles of the distribution implied
#' by \eqn{\log f}{log f} and treating those quantiles as data in a standard
#' Box-Cox analysis.  In the following we use \code{theta} (\eqn{\theta}) to
#' denote the argument of \eqn{\log f}{log f} on the original scale and
#' \code{phi} (\eqn{\phi}) on the Box-Cox transformed scale.
#'
#' @param logf A pointer to a compiled C++ function returning the log
#'   of the target density \eqn{f}.
#' @param ... further arguments to be passed to \code{logf} and related
#'   functions.
#' @param ep_bc A (positive) numeric scalar. Smallest possible value of
#'   \code{phi} to consider.  Used to avoid negative values of \code{phi}.
#' @param min_phi,max_phi  Numeric scalars.  Smallest and largest values
#'   of \code{phi} at which to evaluate \code{logf}, i.e., the range of values
#'   of \code{phi} over which to evaluate \code{logf}. Any components in
#'   \code{min_phi} that are not positive are set to \code{ep_bc}.
#' @param num A numeric scalar. Number of values at which to evaluate
#'   \code{logf}.
#' @param xdiv A numeric scalar.  Only values of \code{phi} at which the
#'   density \eqn{f} is greater than the (maximum of \eqn{f}) / \code{xdiv} are
#'   used.
#' @param probs A numeric scalar. Probabilities at which to estimate the
#'   quantiles of that will be used as data to find lambda.
#' @param lambda_range A numeric vector of length 2.  Range of lambda over
#'   which to optimise.
#' @param phi_to_theta A pointer to a compiled C++ function returning
#'   (the inverse) of the transformation from \code{theta} to \code{phi} used
#'   to ensure positivity of \code{phi} prior to Box-Cox transformation.  The
#'   argument is \code{phi} and the returned value is \code{theta}.  If
#'   \code{phi_to_theta} is undefined at the input value then the function
#'   should return \code{NA}.
#' @param user_args A list of numeric components providing arguments to
#'   the user-supplied functions \code{phi_to_theta} and \code{log_j}.
#' @param log_j A pointer to a compiled C++ function returning the log of
#'  the Jacobian of the transformation from \code{theta} to \code{phi}, i.e.,
#'  based on derivatives of \eqn{phi} with respect to \eqn{theta}. Takes
#'  \code{theta} as its argument.
#'  If this is not supplied then a constant Jacobian is used.
#' @details The general idea is to estimate quantiles of \eqn{f} corresponding
#'   to a set of equally-spaced probabilities in \code{probs} and to use these
#'   estimated quantiles as data in a standard estimation of the Box-Cox
#'   transformation parameter \code{lambda}.
#'
#'   The density \eqn{f} is first evaluated at \code{num} points equally spaced
#'   over the interval (\code{min_phi}, \code{max_phi}).  The continuous
#'   density \eqn{f} is approximated by attaching trapezium-rule estimates of
#'   probabilities to the midpoints of the intervals between the points.  After
#'   standardizing to account for the fact that \eqn{f} may not be normalized,
#'   (\code{min_phi}, \code{max_phi}) is reset so that values with small
#'   estimated probability (determined by \code{xdiv}) are excluded and the
#'   procedure is repeated on this new range.  Then the required quantiles are
#'   estimated by inferring them from a weighted empirical distribution
#'   function based on treating the midpoints as data and the estimated
#'   probabilities at the midpoints as weights.
#' @return A list containing the following components
#'   \item{lambda}{A numeric scalar.  The value of \code{lambda}.}
#'   \item{gm}{A numeric scalar.  Box-cox scaling parameter, estimated by the
#'     geometric mean of the quantiles used in the optimisation to find the
#'     value of lambda.}
#'   \item{init_psi}{A numeric scalar.  An initial estimate of the mode of the
#'     Box-Cox transformed density}
#'   \item{sd_psi}{A numeric scalar.  Estimates of the marginal standard
#'     deviations of the Box-Cox transformed variables.}
#'  \item{phi_to_theta}{as detailed above (only if \code{phi_to_theta} is
#'    supplied)}
#'  \item{log_j}{as detailed above (only if \code{log_j} is supplied)}
#'  \item{user_args}{as detailed above (only if \code{user_args} is supplied)}
#'
#' @references Box, G. and Cox, D. R. (1964) An Analysis of Transformations.
#'  Journal of the Royal Statistical Society. Series B (Methodological), 26(2),
#'  211-252.
#' @references Andrews, D. F. and Gnanadesikan, R. and Warner, J. L. (1971)
#'  Transformations of Multivariate Data, Biometrics, 27(4).
#' @references Eddelbuettel, D. and Francois, R. (2011). Rcpp: Seamless
#'  R and C++ Integration. \emph{Journal of Statistical Software},
#'  \strong{40}(8), 1-18. \doi{10.18637/jss.v040.i08}
#' @references Eddelbuettel, D. (2013). \emph{Seamless R and C++ Integration
#'  with Rcpp}, Springer, New York. ISBN 978-1-4614-6867-7.
#' @examples
#'
#' # Log-normal density ===================
#'
#' # Note: the default value of max_phi = 10 is OK here but this will not
#' # always be the case.
#'
#' ptr_lnorm <- create_xptr("logdlnorm")
#' mu <- 0
#' sigma <- 1
#' lambda <- find_lambda_one_d_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma)
#' lambda
#' x <- ru_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma, log = TRUE, d = 1,
#'              n = 1000, trans = "BC", lambda = lambda)
#'
#' # Gamma density ===================
#'
#' alpha <- 1
#' # Choose a sensible value of max_phi
#' max_phi <- qgamma(0.999, shape = alpha)
#' # [I appreciate that typically the quantile function won't be available.
#' # In practice the value of lambda chosen is quite insensitive to the choice
#' # of max_phi, provided that max_phi is not far too large or far too small.]
#'
#' ptr_gam <- create_xptr("logdgamma")
#' lambda <- find_lambda_one_d_rcpp(logf = ptr_gam, alpha = alpha,
#'                                  max_phi = max_phi)
#' lambda
#' x <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = 1000, trans = "BC",
#'              lambda = lambda)
#'
#' alpha <- 0.1
#' # NB. for alpha < 1 the gamma(alpha, beta) density is not bounded
#' # So the ratio-of-uniforms emthod can't be used but it may work after a
#' # Box-Cox transformation.
#' # find_lambda_one_d() works much better than find_lambda() here.
#'
#' max_phi <- qgamma(0.999, shape = alpha)
#' lambda <- find_lambda_one_d_rcpp(logf = ptr_gam, alpha = alpha,
#'                                  max_phi = max_phi)
#' lambda
#' x <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = 1000, trans = "BC",
#'              lambda = lambda)
#' \donttest{
#' plot(x)
#' plot(x, ru_scale = TRUE)
#' }
#' @seealso \code{\link{ru_rcpp}} to perform ratio-of-uniforms sampling.
#' @seealso \code{\link{find_lambda_rcpp}} to produce (somewhat) automatically
#'   a list for the argument \code{lambda} of \code{ru} for any value of
#'   \code{d}.
#'
#' @export
find_lambda_one_d_rcpp <- function(logf, ..., ep_bc = 1e-4, min_phi = ep_bc,
                              max_phi = 10, num = 1001L, xdiv = 100,
                              probs = seq(0.01, 0.99, by = 0.01),
                              lambda_range = c(-3, 3), phi_to_theta = NULL,
                              log_j = NULL, user_args = list()) {
  # Check that logf is an external pointer.
  is_pointer <- (class(logf) == "externalptr")
  if (!is_pointer) {
    stop("logf must be an external pointer to a function")
  }
  # Check if phi_to_theta and log_j are supplied then they are external
  # pointers.
  is_pointer <- (class(phi_to_theta) == "externalptr")
  if (!is_pointer & !is.null(phi_to_theta)) {
    stop("phi_to_theta must be an external pointer to a function or NULL")
  }
  is_pointer <- (class(log_j) == "externalptr")
  if (!is_pointer & !is.null(log_j)) {
    stop("log_j must be an external pointer to a function or NULL")
  }
  # Check that max_phi > min_phi in all cases
  if (any(max_phi-min_phi <= 0)) {
    stop("max_phi must be larger than min_phi elementwise.")
  }
  if (any(max_phi <= 0)) {
    stop("all components of max_phi must be positive")
  }
  # Set to ep_bc any non-positive elements in min_phi.
  min_phi <- ifelse(min_phi > 0, min_phi, ep_bc)
  # If phi_to_theta is supplied save them to return them later
  trans_list <- list()
  if (!is.null(phi_to_theta)) {
    trans_list$phi_to_theta <- phi_to_theta
    trans_list$log_j <- log_j
    if (is.null(log_j)) {
      log_j = create_log_jac_xptr("log_none_jac")
      trans_list$log_j <- log_j
    }
  } else {
    phi_to_theta = null_phi_to_theta_xptr("no_trans")
    log_j = create_log_jac_xptr("log_none_jac")
  }
  trans_list$user_args <- user_args
  # Set num equally-spaced values of x in [min_phi, max_phi]
  x <- seq(min_phi, max_phi, len = num)
  # Set parameters for passing to cpp_log_rho_4().
  pars <- list(...)
  # Function to determine how deep a list is, i.e. how many layers
  # of listing it has.
  list_depth <- function(x) {
    ifelse(is.list(x), 1L + max(sapply(x, list_depth)), 0L)
  }
  # Find the depth of pars.
  if (length(pars) > 0) {
    pars_depth <- list_depth(pars)
  } else {
    pars_depth <- 0
  }
  # If the user has supplied a list rather than individual components
  # then remove the extra layer of list and retrieve the original
  # variable names.
  if (pars_depth > 1) {
    par_names <- names(pars)
    pars <- unlist(pars, recursive = FALSE)
    # Remove the user's list name, if they gave it one.
    if (!is.null(par_names)) {
      keep_name <- nchar(par_names) + 2
      names(pars) <- substring(names(pars), keep_name)
    }
  }
  tpars <- list()
  ptpfun <- create_psi_to_phi_xptr("no_trans")
  d <- 1
  logf_args <- list(psi_mode = rep(0, d), rot_mat = diag(d), hscale = 0,
                    logf = logf, pars = pars, tpars = tpars, ptpfun = ptpfun,
                    phi_to_theta = phi_to_theta, log_j = log_j,
                    user_args = user_args)
  # Calculate the density (weights) at these values
  log_w <- do.call(rcpp_apply, c(list(x = matrix(x)), logf_args))
  # Shift log_w so that it has a maximum of 0, to try to avoid underflow.
  log_w <- log_w - max(log_w, na.rm = TRUE)
  # Evaluate the density values.
  w <- exp(log_w)
  # Standardise, so that the weights sum to 1
  w <- w / sum(w)
  n <- length(w)
  # Calculate the mean value on each interval between xs
  wbar <- (w[1:(n-1)] + w[2:n]) / 2
  # Calculate the widths of the intervals
  xdiff <- x[2]-x[1]
  # Find the midpoints of the intervals
  xmid <- (x[1:(n-1)] + x[2:n]) / 2
  # Approximate the probability on each interval (by area of trapezium)
  areas <- wbar * xdiff
  # Which areas are greater than the maximum area / xdiv?
  r <- range(xmid[areas > max(areas) / xdiv])
  # Reset the xs over this new range
  x <- seq(r[1], r[2], len = num)
  # Recalculate the weights and the areas and midpoints
  log_w <- do.call(rcpp_apply, c(list(x = matrix(x)), logf_args))
  log_w <- log_w - max(log_w, na.rm = TRUE)
  w <- exp(log_w)
  w <- w / sum(w)
  n <- length(w)
  wbar <- (w[1:(n-1)] + w[2:n]) / 2
  xdiff <- x[2]-x[1]
  xmid <- (x[1:(n-1)] + x[2:n]) / 2
  areas <- wbar * xdiff
  w <- areas / sum(areas)
  # Estimate the 100*probs% quantiles of the density
  qs <- stats::quantile(wecdf(xmid, w), probs = probs)
  qs <- matrix(qs, ncol = 1)
  # Use the quantiles qs as a sample of data to estimate lambda
  # Set to 1 all the input weights associated with the quantiles
  w_ones <- rep(1, nrow(qs))
  temp <- optim_box_cox(x = qs, w = w_ones, lambda_range = lambda_range)
  # Calculate an estimate of the mode of the transformed density
  fy <- w * (xmid / temp$gm) ^ (1 - temp$lambda)
  mode_y <- which.max(fy)
  temp$init_psi <- box_cox(xmid[mode_y], lambda = temp$lambda, gm = temp$gm)
  # Calculate an estimate of the standard deviation of the transformed density
  low_q <- qs[which.min(abs(probs-stats::pnorm(-1 / 2)))]
  up_q <- qs[which.min(abs(probs-stats::pnorm(1 / 2)))]
  temp$sd_psi <- up_q - low_q
  return(c(temp, trans_list))
}

# ========================= find_lambda_rcpp =========================

#' Selecting the Box-Cox parameter for general d using Rcpp
#'
#' Finds a value of the Box-Cox transformation parameter lambda for which
#' the (positive) random variable with log-density \eqn{\log f}{log f} has a
#' density closer to that of a Gaussian random variable.
#' In the following we use \code{theta} (\eqn{\theta}) to denote the argument
#' of \code{logf} on the original scale and \code{phi} (\eqn{\phi}) on the
#' Box-Cox transformed scale.
#'
#' @param logf A pointer to a compiled C++ function returning the log
#'   of the target density \eqn{f}.
#' @param ... further arguments to be passed to \code{logf} and related
#'   functions.
#' @param d A numeric scalar. Dimension of \eqn{f}.
#' @param n_grid A numeric scalar.  Number of ordinates for each variable in
#'   \code{phi}.  If this is not supplied a default value of
#'   \code{ceiling(2501 ^ (1 / d))} is used.
#' @param ep_bc A (positive) numeric scalar. Smallest possible value of
#'   \code{phi} to consider.  Used to avoid negative values of \code{phi}.
#' @param min_phi,max_phi  Numeric vectors.  Smallest and largest values
#'   of \code{phi} at which to evaluate \code{logf}, i.e., the range of values
#'   of \code{phi} over which to evaluate \code{logf}. Any components in
#'   \code{min_phi} that are not positive are set to \code{ep_bc}.
#' @param which_lam A numeric vector.  Contains the indices of the components
#'   of \code{phi} that ARE to be Box-Cox transformed.
#' @param lambda_range A numeric vector of length 2.  Range of lambda over
#'   which to optimise.
#' @param init_lambda A numeric vector of length 1 or \code{d}.  Initial value
#'   of lambda used in the search for the best lambda.  If \code{init_lambda}
#'   is a scalar then \code{rep(init_lambda, d)} is used.
#' @param phi_to_theta A pointer to a compiled C++ function returning
#'   (the inverse) of the transformation from \code{theta} to \code{phi} used
#'   to ensure positivity of \code{phi} prior to Box-Cox transformation.  The
#'   argument is \code{phi} and the returned value is \code{theta}.  If
#'   \code{phi_to_theta} is undefined at the input value then the function
#'   should return \code{NA}.
#' @param user_args A list of numeric components providing arguments to
#'   the user-supplied functions \code{phi_to_theta} and \code{log_j}.
#' @param log_j A pointer to a compiled C++ function returning the log of
#'  the Jacobian of the transformation from \code{theta} to \code{phi}, i.e.,
#'  based on derivatives of \eqn{phi} with respect to \eqn{theta}. Takes
#'  \code{theta} as its argument.
#' @details The general idea is to evaluate the density \eqn{f} on a
#'  \code{d}-dimensional grid, with \code{n_grid} ordinates for each of the
#'  \code{d} variables.
#'  We treat each combination of the variables in the grid as a data point
#'  and perform an estimation of the Box-Cox transformation parameter
#'  \code{lambda}, in which each data point is weighted by the density
#'  at that point.  The vectors \code{min_phi} and \code{max_phi} define the
#'  limits of the grid and \code{which_lam} can be used to specify that only
#'  certain components of \code{phi} are to be transformed.
#' @return A list containing the following components
#'   \item{lambda}{A numeric vector.  The value of \code{lambda}.}
#'   \item{gm}{A numeric vector.  Box-cox scaling parameter, estimated by the
#'     geometric mean of the values of phi used in the optimisation to find
#'     the value of lambda, weighted by the values of f evaluated at phi.}
#'   \item{init_psi}{A numeric vector.  An initial estimate of the mode of the
#'     Box-Cox transformed density}
#'   \item{sd_psi}{A numeric vector.  Estimates of the marginal standard
#'     deviations of the Box-Cox transformed variables.}
#'  \item{phi_to_theta}{as detailed above (only if \code{phi_to_theta} is
#'    supplied)}
#'  \item{log_j}{as detailed above (only if \code{log_j} is supplied)}
#'  \item{user_args}{as detailed above (only if \code{user_args} is supplied)}
#' @references Box, G. and Cox, D. R. (1964) An Analysis of Transformations.
#'  Journal of the Royal Statistical Society. Series B (Methodological), 26(2),
#'  211-252.
#' @references Andrews, D. F. and Gnanadesikan, R. and Warner, J. L. (1971)
#'  Transformations of Multivariate Data, Biometrics, 27(4).
#' @references Eddelbuettel, D. and Francois, R. (2011). Rcpp: Seamless
#'  R and C++ Integration. \emph{Journal of Statistical Software},
#'  \strong{40}(8), 1-18. \doi{10.18637/jss.v040.i08}
#' @references Eddelbuettel, D. (2013). \emph{Seamless R and C++ Integration
#'  with Rcpp}, Springer, New York. ISBN 978-1-4614-6867-7.
#' @examples
#'
#' # Log-normal density ===================
#' # Note: the default value max_phi = 10 is OK here but this will not always
#' # be the case
#' ptr_lnorm <- create_xptr("logdlnorm")
#' mu <- 0
#' sigma <- 1
#' lambda <- find_lambda_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma)
#' lambda
#' x <- ru_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma, d = 1, n = 1000,
#'              trans = "BC", lambda = lambda)
#'
#' # Gamma density ===================
#' alpha <- 1
#' #  Choose a sensible value of max_phi
#' max_phi <- qgamma(0.999, shape = alpha)
#' # [Of course, typically the quantile function won't be available.  However,
#' # In practice the value of lambda chosen is quite insensitive to the choice
#' # of max_phi, provided that max_phi is not far too large or far too small.]
#'
#' ptr_gam <- create_xptr("logdgamma")
#' lambda <- find_lambda_rcpp(logf = ptr_gam, alpha = alpha, max_phi = max_phi)
#' lambda
#' x <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = 1000, trans = "BC",
#'              lambda = lambda)
#'
#' \donttest{
#' # Generalized Pareto posterior distribution ===================
#'
#' n <- 1000
#' # Sample data from a GP(sigma, xi) distribution
#' gpd_data <- rgpd(m = 100, xi = -0.5, sigma = 1)
#' # Calculate summary statistics for use in the log-likelihood
#' ss <- gpd_sum_stats(gpd_data)
#' # Calculate an initial estimate
#' init <- c(mean(gpd_data), 0)
#'
#' n <- 1000
#' # Sample on original scale, with no rotation ----------------
#' ptr_gp <- create_xptr("loggp")
#' for_ru_rcpp <- c(list(logf = ptr_gp, init = init, d = 2, n = n,
#'                      lower = c(0, -Inf)), ss, rotate = FALSE)
#' x1 <- do.call(ru_rcpp, for_ru_rcpp)
#' plot(x1, xlab = "sigma", ylab = "xi")
#' # Parameter constraint line xi > -sigma/max(data)
#' # [This may not appear if the sample is far from the constraint.]
#' abline(a = 0, b = -1 / ss$xm)
#' summary(x1)
#'
#' # Sample on original scale, with rotation ----------------
#' for_ru_rcpp <- c(list(logf = ptr_gp, init = init, d = 2, n = n,
#'                       lower = c(0, -Inf)), ss)
#' x2 <- do.call(ru_rcpp, for_ru_rcpp)
#' plot(x2, xlab = "sigma", ylab = "xi")
#' abline(a = 0, b = -1 / ss$xm)
#' summary(x2)
#'
#' # Sample on Box-Cox transformed scale ----------------
#'
#' # Find initial estimates for phi = (phi1, phi2),
#' # where phi1 = sigma
#' #   and phi2 = xi + sigma / max(x),
#' # and ranges of phi1 and phi2 over over which to evaluate
#' # the posterior to find a suitable value of lambda.
#' temp <- do.call(gpd_init, ss)
#' min_phi <- pmax(0, temp$init_phi - 2 * temp$se_phi)
#' max_phi <- pmax(0, temp$init_phi + 2 * temp$se_phi)
#'
#' # Set phi_to_theta() that ensures positivity of phi
#' # We use phi1 = sigma and phi2 = xi + sigma / max(data)
#'
#' # Create an external pointer to this C++ function
#' ptr_phi_to_theta_gp <- create_phi_to_theta_xptr("gp")
#' # Note: log_j is set to zero by default inside find_lambda_rcpp()
#' lambda <- find_lambda_rcpp(logf = ptr_gp, ss = ss, d = 2, min_phi = min_phi,
#'                            max_phi = max_phi, user_args = list(xm = ss$xm),
#'                            phi_to_theta = ptr_phi_to_theta_gp)
#' lambda
#'
#' # Sample on Box-Cox transformed, without rotation
#' x3 <- ru_rcpp(logf = ptr_gp, ss = ss, d = 2, n = n, trans = "BC",
#'               lambda = lambda, rotate = FALSE)
#' plot(x3, xlab = "sigma", ylab = "xi")
#' abline(a = 0, b = -1 / ss$xm)
#' summary(x3)
#'
#' # Sample on Box-Cox transformed, with rotation
#' x4 <- ru_rcpp(logf = ptr_gp, ss = ss, d = 2, n = n, trans = "BC",
#'               lambda = lambda)
#' plot(x4, xlab = "sigma", ylab = "xi")
#' abline(a = 0, b = -1 / ss$xm)
#' summary(x4)
#'
#' def_par <- graphics::par(no.readonly = TRUE)
#' par(mfrow = c(2,2), mar = c(4, 4, 1.5, 1))
#' plot(x1, xlab = "sigma", ylab = "xi", ru_scale = TRUE,
#'   main = "mode relocation")
#' plot(x2, xlab = "sigma", ylab = "xi", ru_scale = TRUE,
#'   main = "mode relocation and rotation")
#' plot(x3, xlab = "sigma", ylab = "xi", ru_scale = TRUE,
#'   main = "Box-Cox and mode relocation")
#' plot(x4, xlab = "sigma", ylab = "xi", ru_scale = TRUE,
#'   main = "Box-Cox, mode relocation and rotation")
#' graphics::par(def_par)
#' }
#' @seealso \code{\link{ru_rcpp}} to perform ratio-of-uniforms sampling.
#' @seealso \code{\link{find_lambda_one_d_rcpp}} to produce (somewhat)
#'   automatically a list for the argument \code{lambda} of \code{ru} for the
#'   \code{d} = 1 case.
#' @export
find_lambda_rcpp <- function(logf, ..., d = 1, n_grid = NULL, ep_bc = 1e-4,
                        min_phi = rep(ep_bc, d), max_phi = rep(10, d),
                        which_lam = 1:d, lambda_range = c(-3,3),
                        init_lambda = NULL, phi_to_theta = NULL,
                        log_j = NULL, user_args = list()) {
  # Check that logf is an external pointer.
  is_pointer <- (class(logf) == "externalptr")
  if (!is_pointer) {
    stop("logf must be an external pointer to a function")
  }
  # Check if phi_to_theta and log_j are supplied then they are external
  # pointers.
  is_pointer <- (class(phi_to_theta) == "externalptr")
  if (!is_pointer & !is.null(phi_to_theta)) {
    stop("phi_to_theta must be an external pointer to a function or NULL")
  }
  is_pointer <- (class(log_j) == "externalptr")
  if (!is_pointer & !is.null(log_j)) {
    stop("log_j must be an external pointer to a function or NULL")
  }
  #
  if (!is.null(init_lambda)) {
    if (!is.vector(init_lambda)) {
      stop("init_lambda must be a vector.")
    }
    if (length(init_lambda) == 1) {
      init_lambda <- rep(init_lambda, d)
    }
    if (length(init_lambda) != d) {
      stop("init_lambda must be a vector of length 1 or d.")
    }
  }
  # Check that max_phi > min_phi in all cases
  if (any(max_phi-min_phi <= 0)) {
    stop("max_phi must be larger than min_phi elementwise.")
  }
  # Set to ep_bc any non-positive elements in min_phi and max_phi.
  min_phi <- ifelse(min_phi > 0, min_phi, ep_bc)
  max_phi <- ifelse(max_phi > 0, max_phi, ep_bc)
  # If phi_to_theta is supplied save them to return them later
  trans_list <- list()
  if (!is.null(phi_to_theta)) {
    trans_list$phi_to_theta <- phi_to_theta
    trans_list$log_j <- log_j
    if (is.null(log_j)) {
      log_j = create_log_jac_xptr("log_none_jac")
      trans_list$log_j <- log_j
    }
  } else {
    phi_to_theta = null_phi_to_theta_xptr("no_trans")
    log_j = create_log_jac_xptr("log_none_jac")
  }
  trans_list$user_args <- user_args
  # Evaluate the target density (up to a multiplicative constant) over a
  # grid of values that contains most of the probability.
  #
  # If n_grid has not been specified then set a default value
  if (is.null(n_grid)) {
    n_grid <- n_grid_fn(d)
  }
  # Set the minimum and maximum for each variable
  low_phi <- min_phi
  up_phi <- max_phi
  # Create matrix with each column giving the values of each variable
  phi <- mapply(seq, from = low_phi, to = up_phi, len = n_grid)
  # Make this matrix into a list with an entry for each column
  phi <- lapply(seq_len(ncol(phi)), function(i) phi[, i])
  # Expand into a matrix containing the grid of combinations (one
  # combination in each row).
  phi <- expand.grid(phi)
  # Set parameters for passing to cpp_log_rho_4().
  pars <- list(...)
  # Function to determine how deep a list is, i.e. how many layers
  # of listing it has.
  list_depth <- function(x) {
    ifelse(is.list(x), 1L + max(sapply(x, list_depth)), 0L)
  }
  # Find the depth of pars.
  if (length(pars) > 0) {
    pars_depth <- list_depth(pars)
  } else {
    pars_depth <- 0
  }
  # If the user has supplied a list rather than individual components
  # then remove the extra layer of list and retrieve the original
  # variable names.
  if (pars_depth > 1) {
    par_names <- names(pars)
    pars <- unlist(pars, recursive = FALSE)
    # Remove the user's list name, if they gave it one.
    if (!is.null(par_names)) {
      keep_name <- nchar(par_names) + 2
      names(pars) <- substring(names(pars), keep_name)
    }
  }
  tpars <- list()
  ptpfun <- create_psi_to_phi_xptr("no_trans")
  logf_args <- list(psi_mode = rep(0, d), rot_mat = diag(d), hscale = 0,
                    logf = logf, pars = pars, tpars = tpars, ptpfun = ptpfun,
                    phi_to_theta = phi_to_theta, log_j = log_j,
                    user_args = user_args)
  # Evaluate the target log-density at each combination in the grid.
  phi <- data.matrix(phi)
  log_w <- do.call(rcpp_apply, c(list(x = phi), logf_args))
  # Shift log_w so that it has a maximum of 0, to try to avoid underflow.
  log_w <- log_w - max(log_w, na.rm = TRUE)
  # Evaluate the density values.
  w <- exp(log_w)
  #
  # Seek marginal Box-Cox transformations such that the joint posterior is
  # closer to being bivariate normal
  #
  if (any(phi[, which_lam] <= 0)) {
    stop("Attempt to use Box-Cox transformation on non-positive value(s).")
  }
  res_bc <- optim_box_cox(x = phi, w = w, which_lam = which_lam,
                          lambda_range = lambda_range, start = init_lambda)
  lambda <- res_bc$lambda
  gm <- res_bc$gm
  phi_to_psi <- function(phi)  {
    for (j in 1:ncol(phi)) {
      phi_in <- phi[, j]
      phi[, j] <- box_cox(x = phi_in, lambda = lambda[j], gm = gm[j])
    }
    return(phi)
  }
  temp <- init_psi_calc(phi_to_psi = phi_to_psi, phi = phi, lambda = lambda,
                        gm = gm, w = w, which_lam = which_lam)
  temp <- c(list(lambda = lambda, gm = gm), temp)
  return(c(temp, trans_list))
}