# =========================== plot.ru ===========================
#' Plot diagnostics for an ru object
#'
#' \code{plot} method for class \code{"ru"}. For \code{d = 1} a histogram of
#' the simulated values is plotted with a the density function superimposed.
#' The density is normalized crudely using the trapezium rule. For
#' \code{d = 2} a scatter plot of the simulated values is produced with
#' density contours superimposed. For \code{d > 2} pairwise plots of the
#' simulated values are produced.
#'
#' @param x an object of class \code{"ru"}, a result of a call to \code{ru}.
#' @param y Not used.
#' @param ... Additional arguments passed on to \code{hist}, \code{lines},
#' \code{contour} or \code{points}.
#' @param n A numeric scalar. Only relevant if \code{x$d = 1} or
#' \code{x$d = 2}. The meaning depends on the value of x$d.
#' \itemize{
#' \item For d = 1 : n + 1 is the number of abscissae in the trapezium
#' method used to normalize the density.
#' \item For d = 2 : an n by n regular grid is used to contour the density.
#' }
#' @param prob Numeric vector. Only relevant for \code{d = 2}. The contour
#' lines are drawn such that the respective probabilities that the variable
#' lies within the contour are approximately equal to the values in
#' \code{prob}.
#' @param ru_scale A logical scalar. Should we plot data and density on the
#' scale used in the ratio-of-uniforms algorithm (\code{TRUE}) or on the
#' original scale (\code{FALSE})?
#' @param rows A numeric scalar. When \code{d > 2} this sets the number of
#' rows of plots. If the user doesn't provide this then it is set
#' internally.
#' @param xlabs,ylabs Numeric vectors. When \code{d > 2} these set the labels
#' on the x and y axes respectively. If the user doesn't provide these then
#' the column names of the simulated data matrix to be plotted are used.
#' @param var_names A character (or numeric) vector of length \code{x$d}. This
#' argument can be used to replace variable names set using \code{var_names}
#' in the call to \code{\link{ru}} or \code{\link{ru_rcpp}}.
#' @param points_par A list of arguments to pass to
#' \code{\link[graphics]{points}} to control the appearance of points
#' depicting the simulated values. Only relevant when \code{d = 2}.
#' @return No return value, only the plot is produced.
#' @examples
#' # Log-normal density ----------------
#' x <- ru(logf = dlnorm, log = TRUE, d = 1, n = 1000, lower = 0, init = 1)
#' \donttest{
#' plot(x)
#' }
#' # Improve appearance using arguments to plot() and hist()
#' \donttest{
#' plot(x, breaks = seq(0, ceiling(max(x$sim_vals)), by = 0.25),
#' xlim = c(0, 10))
#' }
#' # 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)
#' }
#' x <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1000, init = c(0, 0))
#' \donttest{
#' plot(x)
#' }
#' @seealso \code{\link{summary.ru}} for summaries of the simulated values
#' and properties of the ratio-of-uniforms algorithm.
#' @export
plot.ru <- function(x, y, ..., n = ifelse(x$d == 1, 1001, 101),
prob = c(0.1, 0.25, 0.5, 0.75, 0.95, 0.99),
ru_scale = FALSE, rows = NULL, xlabs = NULL,
ylabs = NULL, var_names = NULL,
points_par = list(col = 8)) {
if (!inherits(x, "ru")) {
stop("use only with \"ru\" objects")
}
if (n < 1) {
stop("n must be no smaller than 1")
}
# Check var_names
if (!is.null(var_names)) {
if (length(var_names) != x$d) {
stop("''var_names'' must have length ''x$d''")
}
colnames(x$sim_vals) <- var_names
}
if (ru_scale) {
plot_data <- x$sim_vals_rho
plot_density <- x$logf_rho
} else {
plot_data <- x$sim_vals
pmedian <- apply(x$sim_vals, 2, stats::median)
hshift <- do.call(x$logf, c(list(pmedian), x$logf_args))
xlogf <- x$logf
plot_density <- function(x, ...) {
xlogf(x, ...) - hshift
}
}
if (!is.null(x$logf_rho_args) & ru_scale) {
density_args <- x$logf_rho_args
} else {
density_args <- x$logf_args
}
if (x$d == 1) {
temp <- graphics::hist(plot_data, plot = FALSE)
a <- temp$breaks[1]
b <- temp$breaks[length(temp$breaks)]
h <- (b - a) / n
xx <- seq(a, b, by = h)
density_fun <- function(z) {
density_list <- c(list(z), density_args)
exp(do.call(plot_density, density_list))
}
yy <- sapply(xx, density_fun)
# Remove any infinite, missing, or undefined values
cond <- is.finite(yy)
yy <- yy[cond]
xx <- xx[cond]
n <- length(yy) - 1
#
area <- h * (yy[1] / 2 + sum(yy[2:n]) + yy[n + 1] / 2)
yy <- yy / area
max_y <- max(temp$density, yy)
temp <- list(...)
my_hist <- function(x, ..., type, lty, lwd, pch, lend, ljoin, lmitre) {
temp_hist <- list(...)
if (is.null(temp_hist$main)) {
graphics::hist(x, main = "", ...)
} else {
graphics::hist(x, ...)
}
}
if (is.null(temp$xlab)) {
my_hist(plot_data, prob = TRUE, ylim = c(0, max_y), xlab = "", ...)
if (!is.null(colnames(plot_data))) {
graphics::title(xlab = parse(text = colnames(plot_data)[1]))
}
} else {
my_hist(plot_data, prob = TRUE, ylim = c(0, max_y), ...)
}
my_lines <- function(x, y, ..., breaks, freq, probability, include.lowest,
right, density, angle, border, plot, labels, nclass) {
graphics::lines(x, y, ...)
}
my_lines(xx, yy, ...)
}
if (x$d == 2) {
r <- apply(plot_data, 2, range)
xx <- seq(r[1, 1], r[2, 1], len = n)
yy <- seq(r[1, 2], r[2, 2], len = n)
zz <- matrix(NA, ncol = length(xx), nrow = length(yy))
for (i in 1:length(xx)) {
for (j in 1:length(yy)) {
for_logf <- c(list(c(xx[i], yy[j])), density_args)
zz[i, j] <- exp(do.call(plot_density, for_logf))
}
}
zz[zz == -Inf] <- NA
dx <- diff(xx[1:2])
dy <- diff(yy[1:2])
sz <- sort(zz)
c1 <- cumsum(sz) * dx * dy
c1 <- c1 / max(c1)
con_levs <- suppressWarnings(sapply(prob, function(x)
stats::approx(c1, sz, xout = 1 - x)$y))
#
graphics::contour(xx, yy, zz, levels = con_levs, add = FALSE, ann = FALSE,
labels = prob * 100, ...)
do.call(graphics::points, c(list(x = plot_data), points_par))
graphics::contour(xx, yy, zz, levels = con_levs, add = TRUE, ann = TRUE,
labels = prob * 100, ...)
temp <- list(...)
if (is.null(temp$xlab)) {
if (!is.null(colnames(plot_data))) {
graphics::title(xlab = parse(text = colnames(plot_data)[1]))
}
}
if (is.null(temp$ylab)) {
if (!is.null(colnames(plot_data))) {
graphics::title(ylab = parse(text = colnames(plot_data)[2]))
}
}
}
if (x$d > 2) {
if (is.null(rows)) {
rows <- x$d - 2
}
cols <- ceiling(choose(x$d, 2) / rows)
if (!is.null(colnames(plot_data))) {
xlabs <- colnames(plot_data)
ylabs <- colnames(plot_data)
} else {
xlabs <- rep(NA, x$d)
ylabs <- rep(NA, x$d)
}
oldpar <- graphics::par(mfrow = c(rows, cols))
on.exit(graphics::par(oldpar))
pairwise_plots <- function(x) {
for (i in 1:(ncol(x) - 1)) {
for (j in (i + 1):ncol(x)) {
graphics::plot(x[, i], x[, j], xlab = "", ylab = "", ...)
graphics::title(xlab = parse(text = xlabs[i]),
ylab = parse(text = ylabs[j]))
}
}
}
pairwise_plots(plot_data)
}
return(invisible())
}
# =========================== summary.ru ===========================
#' Summarizing ratio-of-uniforms samples
#'
#' \code{summary} method for class \code{"ru"}.
#'
#' @param object an object of class \code{"ru"}, a result of a call to
#' \code{ru}.
#' @param ... For \code{summary.lm}: additional arguments passed to
#' \code{\link{summary}}.
#' For \code{print.lm}: additional optional arguments passed to
#' \code{\link{print}}.
#' @return For \code{summary.lm}: a list of the following components from
#' \code{object}:
#' \itemize{
#' \item information about the ratio-of-uniforms bounding box, i.e.,
#' \code{object$box}
#' \item an estimate of the probability of acceptance, i.e.,
#' \code{object$pa}
#' \item a summary of the simulated values, via
#' \code{summary(object$sim_vals)}
#' }
#' @examples
#' # one-dimensional standard normal ----------------
#' x <- ru(logf = function(x) -x ^ 2 / 2, d = 1, n = 1000, init = 0)
#' summary(x)
#'
#' # 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)
#' }
#' x <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1000, init = c(0, 0))
#' summary(x)
#' @seealso \code{\link{ru}} for descriptions of \code{object$sim_vals} and
#' \code{object$box}.
#' @seealso \code{\link{plot.ru}} for a diagnostic plot.
#' @name summary.ru
#' @export
summary.ru <- function(object, ...) {
if (!inherits(object, "ru")) {
stop("use only with \"ru\" objects")
}
sim_summary <- summary(object$sim_vals, ...)
object <- object[c("box", "pa")]
object$sim_summary <- sim_summary
class(object) <- "summary.ru"
return(object)
}
# =========================== print.summary.ru ===========================
#' Print method for objects of class \code{"summary.ru"}
#'
#' \code{print} method for an object \code{object} of class
#' \code{"summary.ru"}.
#'
#' @param x an object of class \code{"summary.ru"}, a result of a call to
#' \code{\link{summary.ru}}.
#' @return For \code{print.summary.ru}: the argument \code{x}, invisibly.
#' @rdname summary.ru
#' @export
print.summary.ru <- function(x, ...) {
if (!inherits(x, "summary.ru")) {
stop("use only with \"summary.ru\" objects")
}
cat("ru bounding box: ", "\n")
print(x$box, ...)
cat("\n")
cat("estimated probability of acceptance: ", "\n")
print(x$pa, ...)
cat("\n")
cat("sample summary", "\n")
print(x$sim_summary, ...)
invisible(x)
}
# ============================= print.ru() ================================== #
#' Print method for an \code{"ru"} object
#'
#' \code{print} method for class \code{"ru"}.
#'
#' @param x an object of class \code{"ru"}, a result of a call to
#' \code{\link{ru}} or \code{\link{ru_rcpp}}.
#' @param ... Additional arguments. None are used in this function.
#' @details Simply prints the call to \code{ru} or \code{ru_rcpp}.
#' @return The argument \code{x}, invisibly.
#' @seealso \code{\link{summary.ru}} for summaries of the simulated values
#' and properties of the ratio-of-uniforms algorithm.
#' @seealso \code{\link{plot.ru}} for a diagnostic plot.
#' @export
print.ru <- function(x, ...) {
if (!inherits(x, "ru")) {
stop("use only with \"ru\" objects")
}
cat("\nCall:\n", paste(deparse(x$call), sep = "\n", collapse = "\n"),
"\n\n", sep = "")
return(invisible(x))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.