R: Generalized ratio-of-uniforms sampling using C++ via Rcpp

ru_rcpp {rust}

R Documentation

Generalized ratio-of-uniforms sampling using C++ via Rcpp

Description

Uses the generalized ratio-of-uniforms method to simulate from a distribution with log-density \log f (up to an additive constant). The density f must be bounded, perhaps after a transformation of variable. The file user_fns.cpp that is sourced before running the examples below is available at the rust Github page at https://raw.githubusercontent.com/paulnorthrop/rust/master/src/user_fns.cpp.

Usage

ru_rcpp(
  logf,
  ...,
  n = 1,
  d = 1,
  init = NULL,
  mode = NULL,
  trans = c("none", "BC", "user"),
  phi_to_theta = NULL,
  log_j = NULL,
  user_args = list(),
  lambda = rep(1L, d),
  lambda_tol = 1e-06,
  gm = NULL,
  rotate = ifelse(d == 1, FALSE, TRUE),
  lower = rep(-Inf, d),
  upper = rep(Inf, d),
  r = 1/2,
  ep = 0L,
  a_algor = if (d == 1) "nlminb" else "optim",
  b_algor = c("nlminb", "optim"),
  a_method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent"),
  b_method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN", "Brent"),
  a_control = list(),
  b_control = list(),
  var_names = NULL,
  shoof = 0.2
)

Arguments

`logf`	An external pointer to a compiled C++ function returning the log of the target density `f` evaluated at its first argument. This function should return `-Inf` when the density is zero. It is better to use `logf =` explicitly, for example, `ru(logf = dnorm, log = TRUE, init = 0.1)`, to avoid argument matching problems. In contrast, `ru(dnorm, log = TRUE, init = 0.1)` will throw an error because partial matching results in `logf` being matched to `log = TRUE`. See the Passing user-supplied C++ functions in the Rcpp Gallery and the Providing a C++ function to `ru_rcpp` section in the Rusting faster: Simulation using Rcpp vignette.
`...`	Further arguments to be passed to `logf` and related functions.
`n`	A non-negative integer scalar. The number of simulated values required. If `n = 0` then no simulation is performed but the component `box` in the returned object gives the ratio-of-uniforms bounding box that would have been used.
`d`	A positive integer scalar. The dimension of `f`.
`init`	A numeric vector. Initial estimates of the mode of `logf`. If `trans = "BC"` or `trans = "user"` this is after Box-Cox transformation or user-defined transformation, but before any rotation of axes. If `init` is not supplied then `rep(1, d)` is used.
`mode`	A numeric vector of length `d`. The mode of `logf`. If `trans = "BC"` or `trans = "user"` this is after Box-Cox transformation or user-defined transformation, but before any rotation of axes. Only supply `mode` if the mode is known: it will not be checked. If `mode` is supplied then `init` is ignored.
`trans`	A character scalar. `trans = "none"` for no transformation, `trans = "BC"` for Box-Cox transformation, `trans = "user"` for a user-defined transformation. If `trans = "user"` then the transformation should be specified using `phi_to_theta` and `log_j` and `user_args` may be used to pass arguments to `phi_to_theta` and `log_j`. See Details and the Examples.
`phi_to_theta`	An external pointer to a compiled C++ function returning (the inverse) of the transformation from `theta` (`\theta`) to `phi` (`\phi`) that may be used to ensure positivity of `\phi` prior to Box-Cox transformation. The argument is `phi` and the returned value is `theta`. If `phi_to_theta` is undefined at the input value then the function should return `NA`. See Details. If `lambda$phi_to_theta` (see argument `lambda` below) is supplied then this is used instead of any function supplied via `phi_to_theta`.
`log_j`	An external pointer to a compiled C++ function returning the log of the Jacobian of the transformation from `theta` (`\theta`) to `phi` (`\phi`), i.e., based on derivatives of `\phi` with respect to `\theta`. Takes `theta` as its argument. If `lambda$log_j` (see argument `lambda` below) is supplied then this is used instead of any function supplied via `log_j`.
`user_args`	A list of numeric components. If trans = ``user'' then `user_args` is a list providing arguments to the user-supplied functions `phi_to_theta` and `log_j`.
`lambda`	Either A numeric vector. Box-Cox transformation parameters, or A list with components lambda A numeric vector. Box-Cox parameters (required). gm A numeric vector. Box-cox scaling parameters (optional). If supplied this overrides any `gm` supplied by the individual `gm` argument described below. init_psi A numeric vector. Initial estimate of mode after Box-Cox transformation (optional). sd_psi A numeric vector. Estimates of the marginal standard deviations of the Box-Cox transformed variables (optional). phi_to_theta as above (optional). log_j As above (optional). user_args As above (optional). This list may be created using `find_lambda_one_d_rcpp` (for `d` = 1) or `find_lambda_rcpp` (for any `d`).
`lambda_tol`	A numeric scalar. Any values in lambda that are less than `lambda_tol` in magnitude are set to zero.
`gm`	A numeric vector. Box-cox scaling parameters (optional). If `lambda$gm` is supplied in input list `lambda` then `lambda$gm` is used, not `gm`.
`rotate`	A logical scalar. If TRUE (`d` > 1 only) use Choleski rotation. If d = 1 and `rotate = TRUE` then rotate will be set to FALSE with a warning. See Details.
`lower`, `upper`	Numeric vectors. Lower/upper bounds on the arguments of the function after any transformation from theta to phi implied by the inverse of `phi_to_theta`. If `rotate = FALSE` these are used in all of the optimisations used to construct the bounding box. If `rotate = TRUE` then they are use only in the first optimisation to maximise the target density.' If `trans = "BC"` components of `lower` that are negative are set to zero without warning and the bounds implied after the Box-Cox transformation are calculated inside `ru`.
`r`	A numeric scalar. Parameter of generalized ratio-of-uniforms.
`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.
`a_algor`, `b_algor`	Character scalars. Either `"nlminb"` or `"optim".` Respective optimisation algorithms used to find `a(r)` and (`b`_i^-(r), `b`_i⁺(r)).
`a_method`, `b_method`	Character scalars. Respective methods used by `optim` to find `a(r)` and (`b`_i^-(r), `b`_i⁺(r)). Only used if `optim` is the chosen algorithm. If `d` = 1 then `a_method` and `b_method` are set to `"Brent"` without warning.
`a_control`, `b_control`	Lists of control arguments to `optim` or `nlminb` to find `a(r)` and (`b`_i^-(r), `b`_i⁺(r)) respectively.
`var_names`	A character (or numeric) vector of length `d`. Names to give to the column(s) of the simulated values.
`shoof`	A numeric scalar in [0, 1]. Sometimes a spurious non-zero convergence indicator is returned from `optim` or `nlminb`). In this event we try to check that a minimum has indeed been found using different algorithm. `shoof` controls the starting value provided to this algorithm. If `shoof = 0` then we start from the current solution. If `shoof = 1` then we start from the initial estimate provided to the previous minimisation. Otherwise, `shoof` interpolates between these two extremes, with a value close to zero giving a starting value that is close to the current solution. The exception to this is when the initial and current solutions are equal. Then we start from the current solution multiplied by `1 - shoof`.

Details

For information about the generalised ratio-of-uniforms method and transformations see the Introducing rust vignette. See also Rusting faster: Simulation using Rcpp These vignettes can also be accessed using vignette("rust-a-vignette", package = "rust") and vignette("rust-c-using-rcpp-vignette", package = "rust").

If trans = "none" and rotate = FALSE then ru implements the (multivariate) generalized ratio of uniforms method described in Wakefield, Gelfand and Smith (1991) using a target density whose mode is relocated to the origin (‘mode relocation’) in the hope of increasing efficiency.

If trans = "BC" then marginal Box-Cox transformations of each of the d variables is performed, with parameters supplied in lambda. The function phi_to_theta may be used, if necessary, to ensure positivity of the variables prior to Box-Cox transformation.

If trans = "user" then the function phi_to_theta enables the user to specify their own transformation.

In all cases the mode of the target function is relocated to the origin after any user-supplied transformation and/or Box-Cox transformation.

If d is greater than one and rotate = TRUE then a rotation of the variable axes is performed after mode relocation. The rotation is based on the Choleski decomposition (see chol) of the estimated Hessian (computed using optimHess of the negated log-density after any user-supplied transformation or Box-Cox transformation. If any of the eigenvalues of the estimated Hessian are non-positive (which may indicate that the estimated mode of logf is close to a variable boundary) then rotate is set to FALSE with a warning. A warning is also given if this happens when d = 1.

The default value of the tuning parameter r is 1/2, which is likely to be close to optimal in many cases, particularly if trans = "BC".

Value

An object of class "ru" is a list containing the following components:

`sim_vals`	An `n` by `d` matrix of simulated values.
`box`	A (2 * `d` + 1) by `d` + 2 matrix of ratio-of-uniforms bounding box information, with row names indicating the box parameter. The columns contain column 1 values of box parameters. columns 2 to (2+`d`-1) values of variables at which these box parameters are obtained. column 2+`d` convergence indicators. Scaling of f within `ru` and relocation of the mode to the origin means that the first row of `box` will always be `c(1, rep(0, d))`.
`pa`	A numeric scalar. An estimate of the probability of acceptance.
`r`	The value of `r`.
`d`	The value of `d`.
`logf`	A function. `logf` supplied by the user, but with f scaled by the maximum of the target density used in the ratio-of-uniforms method (i.e. `logf_rho`), to avoid numerical problems in contouring f in `plot.ru` when `d = 2`.
`logf_rho`	A function. The target function actually used in the ratio-of-uniforms algorithm.
`sim_vals_rho`	An `n` by `d` matrix of values simulated from the function used in the ratio-of-uniforms algorithm.
`logf_args`	A list of further arguments to `logf`.
`logf_rho_args`	A list of further arguments to `logf_rho`. Note: this component is returned by `ru_rcpp` but not by `ru`.
`f_mode`	The estimated mode of the target density f, after any Box-Cox transformation and/or user supplied transformation, but before mode relocation.

References

Wakefield, J. C., Gelfand, A. E. and Smith, A. F. M. (1991) Efficient generation of random variates via the ratio-of-uniforms method. Statistics and Computing (1991), 1, 129-133. doi:10.1007/BF01889987.

Eddelbuettel, D. and Francois, R. (2011). Rcpp: Seamless R and C++ Integration. Journal of Statistical Software, 40(8), 1-18. doi:10.18637/jss.v040.i08

Eddelbuettel, D. (2013). Seamless R and C++ Integration with Rcpp, Springer, New York. ISBN 978-1-4614-6867-7.

Examples

n <- 1000

# Normal density ===================

# One-dimensional standard normal ----------------
ptr_N01 <- create_xptr("logdN01")
x <- ru_rcpp(logf = ptr_N01, d = 1, n = n, init = 0.1)

# Two-dimensional standard normal ----------------
ptr_bvn <- create_xptr("logdnorm2")
rho <- 0
x <- ru_rcpp(logf = ptr_bvn, rho = rho, d = 2, n = n,
  init = c(0, 0))

# Two-dimensional normal with positive association ===================
rho <- 0.9
# No rotation.
x <- ru_rcpp(logf = ptr_bvn, rho = rho, d = 2, n = n, init = c(0, 0),
             rotate = FALSE)

# With rotation.
x <- ru_rcpp(logf = ptr_bvn, rho = rho, d = 2, n = n, init = c(0, 0))

# Using general multivariate normal function.
ptr_mvn <- create_xptr("logdmvnorm")
covmat <- matrix(rho, 2, 2) + diag(1 - rho, 2)
x <- ru_rcpp(logf = ptr_mvn, sigma = covmat, d = 2, n = n, init = c(0, 0))

# Three-dimensional normal with positive association ----------------
covmat <- matrix(rho, 3, 3) + diag(1 - rho, 3)

# No rotation.
x <- ru_rcpp(logf = ptr_mvn, sigma = covmat, d = 3, n = n,
             init = c(0, 0, 0), rotate = FALSE)

# With rotation.
x <- ru_rcpp(logf = ptr_mvn, sigma = covmat, d = 3, n = n,
             init = c(0, 0, 0))

# Log-normal density ===================

ptr_lnorm <- create_xptr("logdlnorm")
mu <- 0
sigma <- 1
# Sampling on original scale ----------------
x <- ru_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma, d = 1, n = n,
             lower = 0, init = exp(mu))

# Box-Cox transform with lambda = 0 ----------------
lambda <- 0
x <- ru_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma, d = 1, n = n,
             lower = 0, init = exp(mu), trans = "BC", lambda = lambda)

# Equivalently, we could use trans = "user" and supply the (inverse) Box-Cox
# transformation and the log-Jacobian by hand
ptr_phi_to_theta_lnorm <- create_phi_to_theta_xptr("exponential")
ptr_log_j_lnorm <- create_log_j_xptr("neglog")
x <- ru_rcpp(logf = ptr_lnorm, mu = mu, sigma = sigma, d = 1, n = n,
  init = 0.1, trans = "user", phi_to_theta = ptr_phi_to_theta_lnorm,
  log_j = ptr_log_j_lnorm)

# Gamma (alpha, 1) density ===================

# Note: the gamma density in unbounded when its shape parameter is < 1.
# Therefore, we can only use trans="none" if the shape parameter is >= 1.

# Sampling on original scale ----------------

ptr_gam <- create_xptr("logdgamma")
alpha <- 10
x <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = n,
  lower = 0, init = alpha)

alpha <- 1
x <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = n,
  lower = 0, init = alpha)

# Box-Cox transform with lambda = 1/3 works well for shape >= 1. -----------

alpha <- 1
x <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = n,
  trans = "BC", lambda = 1/3, init = alpha)
summary(x)

# Equivalently, we could use trans = "user" and supply the (inverse) Box-Cox
# transformation and the log-Jacobian by hand

lambda <- 1/3
ptr_phi_to_theta_bc <- create_phi_to_theta_xptr("bc")
ptr_log_j_bc <- create_log_j_xptr("bc")
x <- ru_rcpp(logf = ptr_gam, alpha = alpha, d = 1, n = n,
  trans = "user", phi_to_theta = ptr_phi_to_theta_bc, log_j = ptr_log_j_bc,
  user_args = list(lambda = lambda), init = alpha)
summary(x)


# Generalized Pareto posterior distribution ===================

# 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
# Mode relocation only ----------------
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)

# Rotation of axes plus mode relocation ----------------
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)

# Cauchy ========================

ptr_c <- create_xptr("logcauchy")

# The bounding box cannot be constructed if r < 1.  For r = 1 the
# bounding box parameters b1-(r) and b1+(r) are attained in the limits
# as x decreases/increases to infinity respectively.  This is fine in
# theory but using r > 1 avoids this problem and the largest probability
# of acceptance is obtained for r approximately equal to 1.26.

res <- ru_rcpp(logf = ptr_c, log = TRUE, init = 0, r = 1.26, n = 1000)

# Half-Cauchy ===================

ptr_hc <- create_xptr("loghalfcauchy")

# Like the Cauchy case the bounding box cannot be constructed if r < 1.
# We could use r > 1 but the mode is on the edge of the support of the
# density so as an alternative we use a log transformation.

x <- ru_rcpp(logf = ptr_hc, init = 0, trans = "BC", lambda = 0, n = 1000)
x$pa
plot(x, ru_scale = TRUE)

# Example 4 from Wakefield et al. (1991) ===================
# Bivariate normal x bivariate student-t

ptr_normt <- create_xptr("lognormt")
rho <- 0.9
covmat <- matrix(c(1, rho, rho, 1), 2, 2)
y <- c(0, 0)

# Case in the top right corner of Table 3
x <- ru_rcpp(logf = ptr_normt, mean = y, sigma1 = covmat, sigma2 = covmat,
  d = 2, n = 10000, init = y, rotate = FALSE)
x$pa

# Rotation increases the probability of acceptance
x <- ru_rcpp(logf = ptr_normt, mean = y, sigma1 = covmat, sigma2 = covmat,
  d = 2, n = 10000, init = y, rotate = TRUE)
x$pa

[Package rust version 1.4.2 Index]