Uses the generalized ratio-of-uniforms (RU) method to simulate
from univariate and (low-dimensional) multivariate continuous distributions.
The user specifies the log-density, up to an additive constant. The RU
algorithm is applied after relocation of mode of the density to zero, and
the user can choose a tuning parameter r. For details see Wakefield, Gelfand
and Smith (1991)
The rust
package implements the multivariate generalized ratio-of-uniforms method of simulating random variates from a d-dimensional continuous distribution. The user specifies (the log of) a positive target function f
that is proportional to the density function of the distribution.
We use the function ru
to simulate a sample of size 1000 from a two-dimensional standard normal distribution with strong positive correlation between the components. Of course, this particular example is purely illustrative: there are better ways to simulate from a multivariate normal distribution.
rho <- 0.9covmat <- 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)}x <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1000, init = c(0, 0))
From version 1.2.0 onwards the faster function ru_rcpp
can be used. See the vignette "Rusting Faster: Simulation using Rcpp" for details.
# Create an external pointer to a C++ function to evaluate the log-density.ptr_bvn <- create_xptr("logdnorm2")# Pass the external pointer to `ru_rcpp`.x <- ru_rcpp(logf = ptr_bvn, rho = rho, d = 2, n = 1000, init = c(0, 0))
To get the current released version from CRAN:
install.packages("rust")
See vignette("rust-vignette", package = "rust")
for an overview of the package, vignette("rust-when-to-use-vignette", package = "rust")
for guidance on when rust
can be used and vignette("rust-using-rcpp-vignette", package = "rust")
for information on how to take advantage of the Rcpp package.
Further steps are taken in the opimisations used to set up the ratio-of-uniforms bounding box (specifically in find_a() and in cpp_find_a()) to avoid erroneous (lack of) convergence warnings. These can occur if we start the optimisation algorithm too close to the solution.
LF line endings used in inst/include/rust.h and inst/include/rust_RcppExports.h to avoid CRAN NOTE
The summary method for class "ru" is now set up according to Section 8.1 of the R FAQ at (https://cran.r-project.org/doc/FAQ/R-FAQ.html).
In the extra checks used to try to avoid erroneous convergence warnings lower and upper bounds on the parameters are now used.
Corrected minor typo in the penultimate sentence of the "Cauchy density" subsection in the "When can rust be used vignette": b1* should be b1+.
In ru()
and ru_rcpp()
the constant hscale
is used to shift logf
(and hence scale the target density f) in function logf
in the returned object. This helps to avoid over/under-flow when contouring f in plot.ru
when d = 2
.
The var_names
argument to ru_rcpp
didn't work. This has been corrected.
The arguments lower, upper
to ru
and ru_rcpp
are now used, at least partly even when rotate = TRUE
. See the updated description of lower, upper
in the documentation of ru
and ru_rcpp
.
That the function logf
supplied to ru
or ru_rcpp
should return -Inf
when the density f is zero is stated explicitly in the help files.
ru() did not work when trans = "user"
and d
> 1. This has been corrected.
Extra checks are used to try to avoid erroneous convergence warnings.
Incorrectly formatted link to the Rcpp Gallery webpage corrected in the vignette "Rusting Faster: Simulation using Rcpp".
Extra examples provided for ru
and for ru_rcpp
: (a) Cauchy, (b) half-Cauchy and (c) bivariate normal x bivariate student-t.
Corrected C++ function vecpow
to avoid compilation errors on some platforms.
Unnecessary dependence on packages devtools
and roxygen2
via Suggests is removed.
Minor edit to vignette: provide link directly to example C++ file user_fns.cpp
in src
directory, rather than the the (identical) example_user_fns.cpp
file in the vignettes
directory.
Packages Rcpp (https://CRAN.R-project.org/package=Rcpp) and RcppArmadillo (https://CRAN.R-project.org/package=RcppArmadillo) are used to speed up the computations if the user provides a C++ function to evaluate their target log-density.
New functions: ru_rcpp
, find_lambda_rcpp
and find_lambda_one_d_rcpp
.
New vignette. "Rusting faster: Simulation using Rcpp".
Bug fixed in plot.ru()
: previously plot.ru()
failed when d > 2
and no axis label names were provided.
Bug fixed in plot.ru
: previously, in the d = 2
case, providing the graphical parameter col
produced an error because col = 8
was hard-coded in a call to points
. Now the extra argument points_par
enables the user to provide a list of arguments to points
.
"using pairs()
" removed from the last sentence Description of plot.ru()
because pairs()
is not used when d > 2
, rather a single plot is produced for each pair of variables.
Obsolete function rho_to_theta()
removed from function ru
in ru_sampling.R
.
If the user calls ru
(or ru_rcpp
) with trans = "user"
but doesn't supply phi_to_theta
then an error is returned.
plot.ru
edited to avoid warning message that occurred in the d=1
case when breaks
was supplied as an argument.
The functions rgpd
, gpd_sum_stats
, gpd_init
and gpd_logpost
are now exported.
plot.ru()
can now be used when d > 2
: it produces pairwise plots of the simulated values.
find_lamba()
: argument init_lambda
added to enable the user to supply an initial value for the Box-Cox transformation vector lambda
.
Unnecessary print statement print(dim(phi))
removed from function find_lambda()
.
Unnecessary print statement print(a_algor)
removed from function ru()
.
Correct lambda$init
to lambda$init_psi
in ru()
when extracting Box-Cox information.
Documentation of ru()
updated to include a description of the returned function logf_rho()
and simulated values sim_vals_rho
and to clarify the meaning of the returned value of f_mode
.
ru()
: the expression for the inverse Box-Cox transformation in the case where lambda is exactly 0 has been corrected.
find_lambda()
: carry out calculation of the target on a shifted log scale to avoid underflow.
Set up plot.ru()
so that if the user supplies axis labels then they are used and otherwise the column name(s) of ru_object$sim_vals
are used. Also enable plotmath symbols to be rendered in the axis labels.