Nothing

```
context("Bounding box")
# Checks the values in the ratio-of-uniforms bounding box returned by ru()
# for some examples in which these values can be found explicitly.
# The examples are:
# A: d-dimensional normal density for d = 1, 2, 3.
# B: 1-dimensional log-normal density under log transformation back to the
# normal density
# C: 1-dimensional gamma density, with shape parameter not less than 1
# A. d-dimensional normal density
normal_box <- function(d, sigma = diag(d), rotate = TRUE, r = 1 / 2) {
#
# Calculates bounding box values in the case of a zero-mean unnormalized
# d-dimensional normal density, with arbitrary covariance structure.
# "Unnormalized" means that the function has a maximum of 1 attained
# at the origin, which is what is required for comparabilty with the
# output from the function ru(), which scales the input function to have
# a maximum of 1 and relocates the mode to the origin.
#
# Args:
# d : dimension of target density.
# sigma : the covariance matrix of the normal density.
# rotate : should the rotation transformation used in the function ru()
# be applied?
# r : ratio-of-uniforms tuning parameter.
#
# Returns:
# box : a (2 * d + 1) by d + 2 matrix of ratio-of-uniforms bounding box
# information with the same structure as object$box returned by
# ru().
#
# Make sure that sigma is a matrix
sigma <- as.matrix(sigma)
# If d > 1 and if rotate = TRUE calculate the covariance matrix after
# the rotation applied in the function ru. The inverse of the
# covariance matrix sigma plays the role of hess_mat in the function ru.
if (rotate & d > 1) {
l_mat <- t(chol(solve(sigma)))
l_mat <- l_mat / det(l_mat) ^ (1/d)
sigma <- t(l_mat) %*% sigma %*% l_mat
}
# Respective values at which bminus and bplus are obtained for each of
# the d marginal variables in the N(0, 1), independent components case.
val <- sqrt((r * d + 1) / r)
# Calculate box
ar <- 1
bplus <- val * exp(-1 /2) * sqrt(diag(sigma))
bminus <- -bplus
box <- c(ar, bminus, bplus)
# Adjustment of val for the marginal variances and for the correlation
# between the components.
adj_mat <- t(sqrt(diag(sigma)) * stats::cov2cor(sigma))
vals <- rbind(rep(0, d), -val * adj_mat, val * adj_mat)
colnames(vals) <- paste("vals", 1:d, sep="")
# We hope to get perfect convegence indicators
conv <- rep(0, 2 * d + 1)
# Create matrix of bounding box information
box <- cbind(box, vals, conv)
bs <- paste(paste("b", 1:d, sep=""),rep(c("minus", "plus"), each=d), sep="")
rownames(box) <- c("a", bs)
#
return(box)
}
my_tol <- 1e-5
# 1. 1-dimensional normal
# (a) N(0, 1)
x1a <- ru(logf = function(x) -x ^ 2 / 2, d = 1, n = 1, init = 0)
test_that("N(0,1)", {
testthat::expect_equal(x1a$box, normal_box(d = 1), tolerance = my_tol)
})
# (b) N(0, 2)
sigma <- 2 # Note: sigma means variance here
x1b <- ru(logf = function(x) -x ^ 2 / (2 * sigma), d = 1, n = 1, init = 0)
test_that("N(0,2)", {
testthat::expect_equal(x1b$box, normal_box(d = 1, sigma = 2),
tolerance = my_tol)
})
# (c) N(1, 2)
sigma <- 2 # Note: sigma means variance here
x1c <- ru(logf = function(x) -(x - 1) ^ 2 / (2 * sigma), d = 1, n = 1,
init = 0)
testthat::expect_equal(x1c$box, normal_box(d = 1, sigma = 2),
tolerance = my_tol)
# 2. 2-dimensional normal
# (a) Zero mean, unit variances and independent components
# Note: rotating shouldn't make any difference in this example.
# (i) rotate = TRUE
x2ai <- ru(logf = function(x) -(x[1]^2 + x[2]^2) / 2, d = 2, n = 1000,
init = c(0, 0), rotate = TRUE)
test_that("BVN, rotation", {
testthat::expect_equal(x2ai$box, normal_box(d = 2), tolerance = my_tol)
})
# (ii) rotate = FALSE
x2aii <- ru(logf = function(x) -(x[1]^2 + x[2]^2) / 2, d = 2, n = 1000,
init = c(0, 0), rotate = FALSE)
test_that("BVN, no rotation", {
testthat::expect_equal(x2aii$box, normal_box(d = 2), tolerance = my_tol)
})
# Function to evaluate the unnormalized log-density of a d-dimensional
# normal distribution
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)
}
# (b) Zero mean, unit variances with positive association
rho <- 0.9
covmat <- matrix(c(1, rho, rho, 1), 2, 2)
# (i) rotate = FALSE
x2bi <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1, init = c(0, 0),
rotate = FALSE)
test_that("BVN, rho = 0.9, no rotation", {
testthat::expect_equal(x2bi$box, normal_box(d = 2, sigma = covmat,
rotate = FALSE),
tolerance = my_tol)
})
# (ii) rotate = TRUE
x2bii <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1, init = c(0, 0),
rotate = TRUE)
test_that("BVN, rho = 0.9, rotation", {
testthat::expect_equal(x2bii$box, normal_box(d = 2, sigma = covmat,
rotate = TRUE),
tolerance = my_tol)
})
# (c) Zero mean, different variances with positive association
covmat <- matrix(c(10, 3, 3, 2), 2, 2)
# (i) rotate = FALSE
x2ci <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1,
init = c(0, 0), rotate = FALSE)
test_that("BVN, general Sigma, no rotation", {
testthat::expect_equal(x2ci$box, normal_box(d = 2, sigma = covmat,
rotate = FALSE),
tolerance = my_tol)
})
# (ii) rotate = TRUE
x2cii <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1,
init = c(0, 0), rotate = TRUE)
test_that("BVN, general Sigma, rotation", {
testthat::expect_equal(x2cii$box, normal_box(d = 2, sigma = covmat,
rotate = TRUE),
tolerance = my_tol)
})
# (d) Mean (1,2), different variances with negative association
covmat <- matrix(c(10, -3, -3, 2), 2, 2)
# (i) rotate = FALSE
x2di <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1,
init = c(0, 0), rotate = FALSE, mean = c(1, 2))
test_that("BVN, non-zero mu, general Sigma, no rotation", {
testthat::expect_equal(x2di$box, normal_box(d = 2, sigma = covmat,
rotate = FALSE),
tolerance = my_tol)
})
# (ii) rotate = TRUE
x2dii <- ru(logf = log_dmvnorm, sigma = covmat, d = 2, n = 1,
init = c(0, 0), rotate = TRUE, mean = c(1, 2))
test_that("BVN, non-zero mu, general Sigma, rotation", {
testthat::expect_equal(x2dii$box, normal_box(d = 2, sigma = covmat,
rotate = TRUE),
tolerance = my_tol)
})
# 3. 3-dimensional normal
# (a) Zero mean, unit variances with positive association
covmat <- matrix(rho, 3, 3) + diag(1 - rho, 3)
# (i) rotate = FALSE
x3ai <- ru(logf = log_dmvnorm, sigma = covmat, d = 3, n = 1,
init = c(0, 0, 0), rotate = FALSE)
test_that("TVN, rho = 0.9, no rotation", {
testthat::expect_equal(x3ai$box, normal_box(d = 3, sigma = covmat,
rotate = FALSE),
tolerance = my_tol)
})
# (ii) rotate = TRUE
x3aii <- ru(logf = log_dmvnorm, sigma = covmat, d = 3, n = 1,
init = c(0, 0, 0), rotate = TRUE)
test_that("TVN, rho = 0.9, rotation", {
testthat::expect_equal(x3aii$box, normal_box(d = 3, sigma = covmat,
rotate = TRUE),
tolerance = my_tol)
})
# B. 1-dimensional log-normal density
# Check that using a log transformation takes us back to the standard
# normal case
lambda <- 0
x <- ru(logf = stats::dlnorm, log = TRUE, d = 1, n = 1, init = 0.1,
trans = "BC", lambda = lambda)
test_that("Log-normal", {
testthat::expect_equal(x$box, normal_box(d = 1), tolerance = my_tol)
})
# C: 1-dimensional gamma density, with shape parameter not less than 1
gamma_box <- function(shape = 1, rate = 1, r = 1 / 2) {
#
# Calculates bounding box values in the case of a zero-mean unnormalized
# d-dimensional normal density, with arbitrary covariance structure.
# "Unnormalized" means that the function has a maximum of 1 attained
# at the origin, which is what is required for comparabilty with the
# output from the function ru(), which scales the input function to have
# a maximum of 1 and relocates the mode to the origin.
#
# Args:
# shape : gamma shape parameter.
# rate : gamma rate parameter.
# r : ratio-of-uniforms tuning parameter.
#
# Returns:
# box : a 3 by 3 matrix of ratio-of-uniforms bounding box information
# with the same structure as object$box returned by ru().
#
if (shape < 1) {
stop("A ratio-of-uniforms bounding box cannot be found when shape < 1")
}
# Mode of the gamma(shape, rate) density
x_mode <- (shape - 1) / rate
# Value of the (normalized) density at the mode
x_max <- stats::dgamma(x_mode, shape = shape, rate = rate)
#
# Calculate box
ar <- 1
a_quad <- rate * r / (r + 1)
b_quad <- -(1 + r * (shape - 1) / (r + 1) - a_quad * x_mode)
c_quad <- -x_mode
b_vals <- Re(polyroot(c(c_quad, b_quad, a_quad)))
vals1 <- c(0, b_vals)
b_fun <- function(x) {
x * (stats::dgamma(x + x_mode, shape = shape) / x_max) ^ (r / (r + 1))
}
box <- c(ar, b_fun(b_vals))
# Create matrix of bounding box information
conv <- rep(0, 3)
box <- cbind(box, vals1, conv)
bs <- paste(paste("b", 1, sep=""),rep(c("minus", "plus"), each=1), sep="")
rownames(box) <- c("a", bs)
#
return(box)
}
# Shape = 1
x1 <- suppressWarnings(ru(logf = dgamma, shape = 1, log = TRUE, d = 1, n = 1,
lower = 0, init = 1))
test_that("Gamma(1, 1)", {
testthat::expect_equal(x1$box, gamma_box(shape = 1), tolerance = my_tol)
})
# Shape = 10
x2 <- ru(logf = dgamma, shape = 10, log = TRUE, d = 1, n = 1, lower = 0,
init = 10)
test_that("Gamma(10, 1)", {
testthat::expect_equal(x2$box, gamma_box(shape = 10), tolerance = my_tol)
})
```

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.