# tests/testthat/test-ru.R In rust: Ratio-of-Uniforms Simulation with Transformation

```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.
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^2 + x^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^2 + x^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)
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)
})
```

## Try the rust package in your browser

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

rust documentation built on June 4, 2021, 1:11 a.m.