Description Usage Arguments Details Value References See Also Examples
View source: R/ru_sampling_rcpp.R
Uses the generalized ratioofuniforms method to simulate from a distribution with logdensity log f (up to an additive constant). f must be bounded, perhaps after a transformation of variable. The file file 'user_fns.cpp' that is sourced before running the examples below is available at the rust Github page at https://github.com/paulnorthrop/rust/blob/master/src/user_fns.cpp.
1 2 3 4 5 6 7 8 9  ru_rcpp(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)

logf 
An external pointer to a compiled C++ 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 
An external pointer to a compiled C++ 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 
An external pointer to a compiled C++ 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 lambda_tol in magnitude are set to zero. 
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. 
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("rustusingrcppvignette", package = "rust")
and
vignette("rustvignette", 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 
logf_rho_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. http://dx.doi.org/10.1007/BF01889987.
Eddelbuettel, D. and Francois, R. (2011). Rcpp: Seamless R and C++ Integration. Journal of Statistical Software, 40(8), 118. http://www.jstatsoft.org/v40/i08/.
Eddelbuettel, D. (2013). Seamless R and C++ Integration with Rcpp, Springer, New York. ISBN 9781461468677.
ru
for a version of ru_rcpp
that
accepts R functions as arguments.
summary.ru
for summaries of the simulated values
and properties of the ratioofuniforms algorithm.
plot.ru
for a diagnostic plot (for d
= 1
and d
= 2 only).
find_lambda_one_d_rcpp
to produce (somewhat)
automatically a list for the argument lambda
of ru
for the
d
= 1 case.
find_lambda_rcpp
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 164 165 166 167 168  n < 1000
# Normal density ===================
# Onedimensional standard normal 
ptr_N01 < create_xptr("logdN01")
x < ru_rcpp(logf = ptr_N01, d = 1, n = n, init = 0.1)
# Twodimensional 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))
# Twodimensional 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))
# Threedimensional 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))
# Lognormal 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))
# BoxCox 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) BoxCox
# transformation and the logJacobian 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)
# BoxCox 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) BoxCox
# transformation and the logJacobian 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)
## Not run:
# 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)
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)
# HalfCauchy ===================
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 studentt
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
## End(Not run)

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