tests/testthat/test-iterative-scheme-mcse.R

In bridgesampling: Bridge Sampling for Marginal Likelihoods and Bayes Factors

context("Iterative scheme MCSE calculations")

test_that("MCSE is finite, positive, and returned for normal method", {
  skip_if_not_installed("Brobdingnag")

  # Minimal deterministic inputs
  q11 <- c(0.10, 0.20, 0.30, 0.40)
  q12 <- c(0.02, 0.09, 0.10, 0.15)
  q21 <- c(0.05, 0.18, 0.22, 0.28)
  q22 <- c(0.01, 0.07, 0.11, 0.12)

  L <- diag(2)
  out <- .run.iterative.scheme(
    q11 = q11,
    q12 = q12,
    q21 = q21,
    q22 = q22,
    r0 = 1,
    tol = 1e-10,
    L = L,
    method = "normal",
    maxiter = 1000,
    silent = TRUE,
    criterion = "r",
    neff = length(q11),
    use_ess = FALSE
  )

  expect_type(out, "list")
  expect_true(is.finite(out$mcse_logml))
  expect_gt(out$mcse_logml, 0)
  # Deterministic value checks for the toy example
  expect_equal(out$logml, 0.13248979618200290, tolerance = 1e-12)
  expect_equal(out$mcse_logml, 0.02353233284705807, tolerance = 1e-12)
})

test_that("MCSE is invariant to constant shifts (warp3 vs normal)", {
  skip_if_not_installed("Brobdingnag")

  # Same base inputs as above
  q11 <- c(0.10, 0.20, 0.30, 0.40)
  q12 <- c(0.02, 0.09, 0.10, 0.15)
  q21 <- c(0.05, 0.18, 0.22, 0.28)
  q22 <- c(0.01, 0.07, 0.11, 0.12)

  L <- diag(2)

  out_normal <- .run.iterative.scheme(
    q11 = q11,
    q12 = q12,
    q21 = q21,
    q22 = q22,
    r0 = 1,
    tol = 1e-10,
    L = L,
    method = "normal",
    maxiter = 1000,
    silent = TRUE,
    criterion = "r",
    neff = length(q11),
    use_ess = FALSE
  )

  out_warp3 <- .run.iterative.scheme(
    q11 = q11,
    q12 = q12,
    q21 = q21,
    q22 = q22,
    r0 = 1,
    tol = 1e-10,
    L = L,
    method = "warp3",
    maxiter = 1000,
    silent = TRUE,
    criterion = "r",
    neff = length(q11),
    use_ess = FALSE
  )

  # warp3 adds a constant to both l1 and l2 and shifts l*, the e^(l - l*) terms are invariant
  expect_equal(out_warp3$mcse_logml, out_normal$mcse_logml, tolerance = 1e-10)
})

test_that("MCSE roughly scales like 1/sqrt(n)", {
  skip_if_not_installed("Brobdingnag")

  base_q11 <- c(0.10, 0.20, 0.30, 0.40)
  base_q12 <- c(0.02, 0.09, 0.10, 0.15)
  base_q21 <- c(0.05, 0.18, 0.22, 0.28)
  base_q22 <- c(0.01, 0.07, 0.11, 0.12)
  L <- diag(2)

  # n
  out_n <- .run.iterative.scheme(
    q11 = base_q11,
    q12 = base_q12,
    q21 = base_q21,
    q22 = base_q22,
    r0 = 1,
    tol = 1e-10,
    L = L,
    method = "normal",
    maxiter = 1000,
    silent = TRUE,
    criterion = "r",
    neff = length(base_q11),
    use_ess = FALSE
  )

  # 4n (replicate samples 4x)
  k <- 4
  out_4n <- .run.iterative.scheme(
    q11 = rep(base_q11, k),
    q12 = rep(base_q12, k),
    q21 = rep(base_q21, k),
    q22 = rep(base_q22, k),
    r0 = 1,
    tol = 1e-10,
    L = L,
    method = "normal",
    maxiter = 1000,
    silent = TRUE,
    criterion = "r",
    neff = length(base_q11) * k,
    use_ess = FALSE
  )

  # Expect MCSE to drop by ~ 1/sqrt(k)
  expect_lt(out_4n$mcse_logml, out_n$mcse_logml)
  expect_equal(out_4n$mcse_logml, out_n$mcse_logml / sqrt(k), tolerance = 0.05)

  # Deterministic value for the 4x replicated case
  expect_equal(out_4n$mcse_logml, 0.01052514482181162, tolerance = 1e-12)
})

test_that("Function runs with use_ess = TRUE (if posterior installed)", {
  skip_if_not_installed("Brobdingnag")
  # if (!requireNamespace("posterior", quietly = TRUE)) {
  #   skip("posterior not installed")
  # } ## uses coda currently

  q11 <- c(0.10, 0.20, 0.30, 0.40, 0.45, 0.5)
  q12 <- c(0.02, 0.09, 0.10, 0.15, 0.18, 0.2)
  q21 <- c(0.05, 0.18, 0.22, 0.28, 0.30, 0.35)
  q22 <- c(0.01, 0.07, 0.11, 0.12, 0.14, 0.15)
  L <- diag(2)

  out <- .run.iterative.scheme(
    q11 = q11,
    q12 = q12,
    q21 = q21,
    q22 = q22,
    r0 = 1,
    tol = 1e-10,
    L = L,
    method = "normal",
    maxiter = 1000,
    silent = TRUE,
    criterion = "r",
    neff = length(q11),
    use_ess = TRUE
  )

  expect_true(is.finite(out$mcse_logml))
  expect_gt(out$mcse_logml, 0)
})