tests/testthat/test-find_lambda.R

In rust: Ratio-of-Uniforms Simulation with Transformation

#context("find_lambda")

# Checks that find_lambda() and find_lambda_rcpp() return the same values.

my_tol <- 1e-5

# 1. Gamma density

alpha <- 1
max_phi <- qgamma(0.999, shape = alpha)
lambda <- find_lambda(logf = dgamma, shape = alpha, log = TRUE,
                      max_phi = max_phi)
lambda
ptr_gam <- create_xptr("logdgamma")
lambda_rcpp <- find_lambda_rcpp(logf = ptr_gam, alpha = alpha,
                                max_phi = max_phi)
lambda_rcpp

testthat::expect_equal(lambda$lambda, lambda_rcpp$lambda, tolerance = my_tol)
testthat::expect_equal(lambda$gm, lambda_rcpp$gm, tolerance = my_tol)
testthat::expect_equal(lambda$init_psi, lambda_rcpp$init_psi,
                       tolerance = my_tol)
testthat::expect_equal(lambda$sd_psi, lambda_rcpp$sd_psi, tolerance = my_tol)

# 2. Posterior density of Generalized Pareto parameters

# Sample data from a GP(sigma, xi) distribution
set.seed(47)
gpd_data <- rgpd(m = 100, xi = 0, sigma = 1)
# Calculate summary statistics for use in the log-likelihood
ss <- gpd_sum_stats(gpd_data)
# Calculate an initial estimate for phi = (phi1, phi2)
temp <- do.call(gpd_init, ss)
min_phi <- pmax(0, temp$init_phi - 2 * temp$se_phi)
max_phi <- pmax(0, temp$init_phi + 2 * temp$se_phi)

# find_lambda() -------------

phi_to_theta <- function(phi) c(phi[1], phi[2] - phi[1] / ss$xm)
log_j <- function(x) 0

lambda <- find_lambda(logf = gpd_logpost, ss = ss, d = 2, min_phi = min_phi,
                      max_phi = max_phi, phi_to_theta = phi_to_theta,
                      log_j = log_j)

# find_lambda_rcpp() -------------

ptr_gp <- create_xptr("loggp")
ptr_phi_to_theta_gp <- create_phi_to_theta_xptr("gp")
log_j <- create_log_jac_xptr("log_none_jac")
lambda_rcpp <- find_lambda_rcpp(logf = ptr_gp, ss = ss, d = 2,
                                min_phi = min_phi, max_phi = max_phi,
                                user_args = list(xm = ss$xm), log_j = log_j,
                                phi_to_theta = ptr_phi_to_theta_gp)

test_that("lambda", {
  testthat::expect_equal(lambda$lambda, lambda_rcpp$lambda, tolerance = my_tol)
})
test_that("gm", {
  testthat::expect_equal(lambda$gm, lambda_rcpp$gm, tolerance = my_tol)
})
test_that("init_psi", {
  testthat::expect_equal(lambda$init_psi, lambda_rcpp$init_psi,
                         tolerance = my_tol)
})
test_that("sd_psi", {
  testthat::expect_equal(lambda$sd_psi, lambda_rcpp$sd_psi, tolerance = my_tol)
})