Description Usage Arguments Details Value References See Also Examples
Uses the generalized ratioofuniforms method to simulate from a distribution with logdensity log f (up to an additive constant). The density f must be bounded, perhaps after a transformation of variable.
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 = 1e06,
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("NelderMead", "BFGS", "CG", "LBFGSB", "SANN", "Brent"),
b_method = c("NelderMead", "BFGS", "CG", "LBFGSB", "SANN", "Brent"),
a_control = list(),
b_control = list(),
var_names = NULL,
shoof = 0.2
)

logf 
A function returning the log of the target density f.
This function should return 
... 
Further arguments to be passed to 
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 
trans 
A character scalar. "none" for no transformation, "BC" for
BoxCox transformation, "user" for a userdefined transformation.
If 
phi_to_theta 
A function returning (the inverse) of the transformation
from theta to phi used to ensure positivity of phi prior to BoxCox
transformation. The argument is phi and the returned value is theta.
If 
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 
lambda 
Either

lambda_tol 
A numeric scalar. Any values in lambda that are less
than 
gm 
A numeric vector. Boxcox scaling parameters (optional). If

rotate 
A logical scalar. If TRUE ( 
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 
r 
A numeric scalar. Parameter of generalized ratioofuniforms. 
ep 
A numeric scalar. Controls initial estimates for optimizations
to find the bbounding box parameters. The default ( 
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

a_control, b_control 
Lists of control arguments to 
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
nonzero convergence indicator is returned from

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 BoxCox 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 BoxCox
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 usersupplied transformation and/or BoxCox 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
logdensity after any usersupplied transformation or BoxCox
transformation. If any of the eigenvalues of the estimated Hessian are
nonpositive (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("rustavignette", package = "rust")
for full details.
An object of class "ru" is a list containing the following components:
sim_vals 
An 
box 
A (2 *
Scaling of f within 
pa 
A numeric scalar. An estimate of the probability of acceptance. 
d 
A numeric scalar. The dimension of 
logf 
A function. 
logf_rho 
A function. The target function actually used in the ratioofuniforms algorithm. 
sim_vals_rho 
An 
logf_args 
A list of further arguments to 
f_mode 
The estimated mode of the target density f, after any BoxCox transformation and/or user supplied transformation, but before mode relocation. 
Wakefield, J. C., Gelfand, A. E. and Smith, A. F. M. (1991) Efficient generation of random variates via the ratioofuniforms method. Statistics and Computing (1991), 1, 129133. 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 ratioofuniforms 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.
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 ===================
# Onedimensional standard normal 
x < ru(logf = function(x) x ^ 2 / 2, d = 1, n = 1000, init = 0.1)
# Twodimensional standard normal 
x < ru(logf = function(x) (x[1]^2 + x[2]^2) / 2, d = 2, n = 1000,
init = c(0, 0))
# Twodimensional 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))
# threedimensional 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))
# Lognormal density ===================
# Sampling on original scale 
x < ru(logf = dlnorm, log = TRUE, d = 1, n = 1000, lower = 0, init = 1)
# BoxCox 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) BoxCox
# transformation and the logJacobian 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)
# BoxCox 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) BoxCox
# transformation and the logJacobian 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 loglikelihood
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)
# HalfCauchy ===================
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 studentt
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

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.