ru: Generalized ratio-of-uniforms sampling In rust: Ratio-of-Uniforms Simulation with Transformation

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.

Usage

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27``` ```ru( logf, ..., n = 1, d = 1, init = 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` A function returning the log of the target density f. This function should return `-Inf` when the density is zero. `...` Further arguments to be passed to `logf` and related functions. `n` A numeric scalar. Number of simulated values required. `d` A numeric scalar. 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. `trans` A character scalar. "none" for no transformation, "BC" for Box-Cox transformation, "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`. `phi_to_theta` A function returning (the inverse) of the transformation from theta to phi 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. `log_j` A function returning the log of the Jacobian of the transformation from theta to phi, i.e. based on derivatives of phi with respect to theta. Takes theta as its argument. `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 lambdaA numeric vector. Box-Cox parameters (required). gmA numeric vector. Box-cox scaling parameters (optional). If supplied this overrides any `gm` supplied by the individual `gm` argument described below. init_psiA numeric vector. Initial estimate of mode after Box-Cox transformation (optional). sd_psiA numeric vector. Estimates of the marginal standard deviations of the Box-Cox transformed variables (optional). phi_to_thetaas above (optional). log_jas above (optional). This list may be created using find_lambda_one_d (for `d` = 1) or find_lambda (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. `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 optimizations 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 optimizations 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 optimization algorithms used to find a(r) and (bi-(r), bi+(r)). `a_method, b_method` Character scalars. Respective methods used by `optim` to find a(r) and (bi-(r), bi+(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 (bi-(r), bi+(r)) respectively. `var_names` A character vector. 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

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"`.

See `vignette("rust-a-vignette", package = "rust")` for full details.

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 1values 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. `d` A numeric scalar. The dimension of `logf`. `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 fin `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`. `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.

`ru_rcpp` for a version of `ru` that uses the Rcpp package to improve efficiency.

`summary.ru` for summaries of the simulated values and properties of the ratio-of-uniforms algorithm.

`plot.ru` for a diagnostic plot.

`find_lambda_one_d` to produce (somewhat) automatically a list for the argument `lambda` of `ru` for the `d` = 1 case.

`find_lambda` to produce (somewhat) automatically a list for the argument `lambda` of `ru` for any value of `d`.

`optim` for choices of the arguments `a_method`, `b_method`, `a_control` and `b_control`.

`nlminb` for choices of the arguments `a_control` and `b_control`.

`optimHess` for Hessian estimation.

`chol` for the Choleski decomposition.

Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163``` ```# Normal density =================== # One-dimensional standard normal ---------------- x <- ru(logf = function(x) -x ^ 2 / 2, d = 1, n = 1000, init = 0.1) # Two-dimensional standard normal ---------------- x <- ru(logf = function(x) -(x[1]^2 + x[2]^2) / 2, d = 2, n = 1000, init = c(0, 0)) # Two-dimensional normal with positive association ---------------- rho <- 0.9 covmat <- matrix(c(1, rho, rho, 1), 2, 2) log_dmvnorm <- function(x, mean = rep(0, d), sigma = diag(d)) { x <- matrix(x, ncol = length(x)) d <- ncol(x) - 0.5 * (x - mean) %*% solve(sigma) %*% t(x - mean) } # No rotation. x <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1000, init = c(0, 0), rotate = FALSE) # With rotation. x <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1000, init = c(0, 0)) # three-dimensional normal with positive association ---------------- covmat <- matrix(rho, 3, 3) + diag(1 - rho, 3) # No rotation. Slow ! x <- ru(logf = log_dmvnorm, sigma = covmat, d = 3, n = 1000, init = c(0, 0, 0), rotate = FALSE) # With rotation. x <- ru(logf = log_dmvnorm, sigma = covmat, d = 3, n = 1000, init = c(0, 0, 0)) # Log-normal density =================== # Sampling on original scale ---------------- x <- ru(logf = dlnorm, log = TRUE, d = 1, n = 1000, lower = 0, init = 1) # Box-Cox transform with lambda = 0 ---------------- lambda <- 0 x <- ru(logf = dlnorm, log = TRUE, d = 1, n = 1000, lower = 0, init = 0.1, trans = "BC", lambda = lambda) # Equivalently, we could use trans = "user" and supply the (inverse) Box-Cox # transformation and the log-Jacobian by hand x <- ru(logf = dlnorm, log = TRUE, d = 1, n = 1000, init = 0.1, trans = "user", phi_to_theta = function(x) exp(x), log_j = function(x) -log(x)) # 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 ---------------- alpha <- 10 x <- ru(logf = dgamma, shape = alpha, log = TRUE, d = 1, n = 1000, lower = 0, init = alpha) alpha <- 1 x <- ru(logf = dgamma, shape = alpha, log = TRUE, d = 1, n = 1000, lower = 0, init = alpha) # Box-Cox transform with lambda = 1/3 works well for shape >= 1. ----------- alpha <- 1 x <- ru(logf = dgamma, shape = alpha, log = TRUE, d = 1, n = 1000, 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 # Note: when phi_to_theta is undefined at x this function returns NA phi_to_theta <- function(x, lambda) { ifelse(x * lambda + 1 > 0, (x * lambda + 1) ^ (1 / lambda), NA) } log_j <- function(x, lambda) (lambda - 1) * log(x) lambda <- 1/3 x <- ru(logf = dgamma, shape = alpha, log = TRUE, d = 1, n = 1000, trans = "user", phi_to_theta = phi_to_theta, log_j = log_j, 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) # Mode relocation only ---------------- n <- 1000 x1 <- ru(logf = gpd_logpost, ss = ss, d = 2, n = n, init = init, lower = c(0, -Inf), rotate = FALSE) 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 ---------------- x2 <- ru(logf = gpd_logpost, ss = ss, d = 2, n = n, init = init, lower = c(0, -Inf)) plot(x2, xlab = "sigma", ylab = "xi") abline(a = 0, b = -1 / ss\$xm) summary(x2) # Cauchy ======================== # 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(logf = dcauchy, log = TRUE, init = 0, r = 1.26, n = 1000) # Half-Cauchy =================== log_halfcauchy <- function(x) { return(ifelse(x < 0, -Inf, dcauchy(x, log = TRUE))) } # 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(logf = log_halfcauchy, 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 log_norm_t <- function(x, mean = rep(0, d), sigma1 = diag(d), sigma2 = diag(d)) { x <- matrix(x, ncol = length(x)) d <- ncol(x) log_h1 <- -0.5 * (x - mean) %*% solve(sigma1) %*% t(x - mean) log_h2 <- -2 * log(1 + 0.5 * x %*% solve(sigma2) %*% t(x)) return(log_h1 + log_h2) } 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(logf = log_norm_t, mean = y, sigma1 = covmat, sigma2 = covmat, d = 2, n = 10000, init = y, rotate = FALSE) x\$pa # Rotation increases the probability of acceptance x <- ru(logf = log_norm_t, mean = y, sigma1 = covmat, sigma2 = covmat, d = 2, n = 10000, init = y, rotate = TRUE) x\$pa ```

rust documentation built on June 4, 2021, 1:11 a.m.