tests/testthat/test-lugsail.R

In LugsailGR: Generalized Gelman-Rubin Diagnostic and Effective Sample Size for MCMC

test_that("calc_ess_bound and calc_psrf_cutoff match theoretical values", {
  # Example 1 from Vats & Knudson (2021): p = 1, alpha = 0.05, epsilon = 0.10, m = 3
  M_val <- calc_ess_bound(p = 1, alpha = 0.05, epsilon = 0.10)
  expect_equal(round(M_val), 1537)

  delta_eps <- calc_psrf_cutoff(m = 3, M_val = M_val)
  expect_equal(round(delta_eps, 6), 1.000976)
})

test_that("lugsail_gr works for univariate and multivariate chains", {
  set.seed(42)
  # Univariate 3 chains
  c1 <- rnorm(500)
  c2 <- rnorm(500)
  c3 <- rnorm(500)
  fit_uni <- lugsail_gr(list(c1, c2, c3))

  expect_s3_class(fit_uni, "lugsail_gr")
  expect_true(fit_uni$psrf[1] > 0)
  expect_true(fit_uni$ess[1] > 0)
  expect_true(is.logical(fit_uni$converged))

  # Multivariate 2 chains
  m1 <- matrix(rnorm(600), ncol = 2)
  m2 <- matrix(rnorm(600), ncol = 2)
  fit_multi <- lugsail_gr(list(m1, m2))

  expect_s3_class(fit_multi, "lugsail_gr")
  expect_true(fit_multi$mpsrf > 0)
  expect_true(fit_multi$mess > 0)
})

test_that("mcmc_gr_test works with user-supplied target PDF and data", {
  set.seed(123)
  obs_data <- rnorm(30, mean = 2, sd = 1)
  log_post <- function(mu, data) {
    sum(dnorm(data, mean = mu, sd = 1, log = TRUE))
  }

  fit <- mcmc_gr_test(target_pdf = log_post, data = obs_data, start_par = c(mu = 0),
                      n_iter = 300, n_chains = 2)

  expect_s3_class(fit, "lugsail_gr")
  expect_true(!is.na(fit$means["mu"]))
})

test_that("censor_gr_test works for right-censored data", {
  set.seed(123)
  times <- rexp(25, rate = 0.5)
  c_time <- 2.0
  obs <- pmin(times, c_time)
  st <- as.numeric(times <= c_time)

  f_exp <- function(x, rate) dexp(x, rate = rate)
  F_exp <- function(x, rate) pexp(x, rate = rate)

  fit <- censor_gr_test(pdf_fn = f_exp, cdf_fn = F_exp, data = obs,
                        start_par = c(rate = 0.6), censor_type = "right",
                        status = st, n_iter = 300, n_chains = 2)

  expect_s3_class(fit, "lugsail_gr")
  expect_true(fit$means["rate"] > 0)
})