tests/testthat/test-validation.R

In AdMit: Adaptive Mixture of Student-t Distributions

## Regressions for the defects found in the 2.1.10 pre-submission review.
## All of them predate 2.1.10.

test_that("fn.optp returns H probabilities when both optimizers fail", {
  ## the fallback assigned the length-H probability vector to the optimizer's
  ## parameter slot, and the logit transform then returned H + 1 of them
  optp <- get("fn.optp", envir = asNamespace("AdMit"))
  for (H in 2:4) {
    set.seed(60 + H)
    Np <- 30L
    lnK <- matrix(rnorm(Np * H), Np, H)
    lnD <- matrix(rnorm(Np * H * H, -1, 0.5), Np, H * H)
    ## iter.max = 1 exhausts both nlminb() and optim(), reaching method "NONE"
    out <- optp(rep(1 / (H - 1), H - 1), lnK, lnD,
                list(trace = 0, iter.max = 1, rel.tol = 1e-8, weightNC = 0.1))
    expect_identical(out$method, "NONE")
    expect_length(out$p, H)
    expect_equal(sum(out$p), 1)
    ## "keep past values": the starting probabilities, unchanged
    expect_equal(out$p, c(0.9 * rep(1 / (H - 1), H - 1), 0.1))
  }
})

test_that("AdMit survives an optimizer that never converges", {
  ## used to abort with "incorrect number of probabilities"
  set.seed(601)
  expect_no_error(AdMit(KERNEL = GelmanMeng, mu0 = c(0.0, 0.1),
                        control = list(Ns = 5e3, Np = 5e2, Hmax = 3, maxit.p = 1)))
})

test_that("dMit(log = TRUE) stays finite where the components are finite", {
  mit <- list(p = c(0.5, 0.5), mu = rbind(c(0, 0), c(1, 1)),
              Sigma = rbind(c(1, 0, 0, 1), c(1, 0, 0, 1)), df = 5)
  theta <- rbind(c(1e10, 1e10), c(1e50, 1e50), c(1e100, 1e100))

  ref <- vapply(seq_len(nrow(theta)), function(i) {
    l <- vapply(seq_along(mit$p), function(h)
      log(mit$p[h]) + mvtnorm::dmvt(theta[i, , drop = FALSE], mit$mu[h, ],
                                    matrix(mit$Sigma[h, ], 2, 2), mit$df, log = TRUE),
      numeric(1))
    max(l) + log(sum(exp(l - max(l))))
  }, numeric(1))

  got <- dMit(theta, mit, log = TRUE)
  expect_true(all(is.finite(got)))
  expect_equal(got, ref)

  ## further out mvtnorm::dmvt itself overflows, so every component really is
  ## -Inf and dMit can do no better than report it
  far <- vapply(seq_along(mit$p), function(h)
    mvtnorm::dmvt(matrix(c(1e200, 1e200), 1), mit$mu[h, ],
                  matrix(mit$Sigma[h, ], 2, 2), mit$df, log = TRUE), numeric(1))
  expect_true(all(!is.finite(far)))
  expect_identical(dMit(c(1e200, 1e200), mit, log = TRUE), -Inf)

  ## univariate, single component, one degree of freedom
  mit1 <- list(p = 1, mu = as.matrix(0), Sigma = as.matrix(1), df = 1)
  expect_true(is.finite(dMit(1e120, mit1, log = TRUE)))
  expect_equal(dMit(1e120, mit1, log = TRUE),
               mvtnorm::dmvt(matrix(1e120), 0, matrix(1), 1, log = TRUE))
})

test_that("dMit leaves representable values exactly as they were", {
  ## the log-scale path is a fall-back, not a rewrite: ordinary values must
  ## keep the result of the plain sum
  mit <- list(p = c(0.3, 0.7), mu = rbind(c(0, 0), c(2, -1)),
              Sigma = rbind(c(1, 0, 0, 1), c(0.5, 0, 0, 0.5)), df = 5)
  set.seed(602)
  theta <- matrix(rnorm(400), 200, 2)
  naive <- 0
  for (h in seq_along(mit$p))
    naive <- naive + exp(log(mit$p[h]) +
                         mvtnorm::dmvt(theta, mit$mu[h, ], matrix(mit$Sigma[h, ], 2, 2),
                                       mit$df, log = TRUE))
  expect_identical(dMit(theta, mit, log = TRUE), log(naive))
})

test_that("AdMitMH works with the documented default mixture", {
  ## AdMitMH read ncol(mit$mu) from the still-empty list and failed with
  ## "invalid 'length' argument"
  K1 <- function(x, log = TRUE) {
    x <- as.matrix(x)
    r <- -0.5 * x[,1]^2
    if (!log) r <- exp(r)
    as.vector(r)
  }
  set.seed(603)
  expect_warning(out <- AdMitMH(N = 500, KERNEL = K1), "not well defined")
  expect_identical(dim(out$draws), c(500L, 1L))
  expect_true(out$accept > 0 && out$accept <= 1)
})

test_that("AdMitIS forwards arguments to a KERNEL or G that takes dots", {
  ## an argument destined for a callee's '...' matched no literal formal and
  ## was dropped, silently changing the estimate
  Kdots <- function(x, log = TRUE, ...) {
    a <- list(...)$aa
    if (is.null(a)) a <- 1
    if (is.vector(x)) x <- matrix(x, nrow = 1)
    r <- -0.5 * a * ((x[,1] - 3)^2 + (x[,2] - 3)^2)
    if (!log) r <- exp(r)
    as.vector(r)
  }
  Kformal <- function(x, aa = 1, log = TRUE) {
    if (is.vector(x)) x <- matrix(x, nrow = 1)
    r <- -0.5 * aa * ((x[,1] - 3)^2 + (x[,2] - 3)^2)
    if (!log) r <- exp(r)
    as.vector(r)
  }
  mit <- list(p = 1, mu = matrix(c(3, 3), 1, 2),
              Sigma = matrix(c(1, 0, 0, 1), 1, 4), df = 5)

  set.seed(605); formal <- AdMitIS(N = 2e4, KERNEL = Kformal, mit = mit, aa = 9)
  set.seed(605); dots <- expect_no_warning(AdMitIS(N = 2e4, KERNEL = Kdots, mit = mit, aa = 9))
  expect_equal(dots$ghat, formal$ghat)
  ## and the argument really did reach the kernel: aa = 1 gives a different answer
  set.seed(605); other <- AdMitIS(N = 2e4, KERNEL = Kformal, mit = mit, aa = 1)
  expect_false(isTRUE(all.equal(dots$ghat, other$ghat)))

  ## the same for G
  Gdots <- function(theta, ...) theta - list(...)$shiftvalue
  set.seed(606)
  shifted <- expect_no_warning(AdMitIS(N = 1e3, KERNEL = Kformal, mit = mit, aa = 9,
                                       G = Gdots, shiftvalue = 1))
  set.seed(606)
  plain <- AdMitIS(N = 1e3, KERNEL = Kformal, mit = mit, aa = 9)
  expect_equal(shifted$ghat, plain$ghat - 1)

  ## and an argument no one can take still warns
  expect_warning(AdMitIS(N = 1e3, KERNEL = Kformal, mit = mit, notanarg = 1),
                 "used by neither")
})

test_that("a constant function of interest gives RNE = NA, not NaN", {
  set.seed(607)
  fit <- AdMit(KERNEL = GelmanMeng, mu0 = c(0.0, 0.1),
               control = list(Ns = 5e3, Np = 5e2, Hmax = 2))
  zero <- AdMitIS(N = 1e3, KERNEL = GelmanMeng, mit = fit$mit,
                  G = function(theta) rep(0, nrow(theta)))
  expect_equal(zero$ghat, 0)
  expect_identical(zero$NSE, 0)
  expect_identical(zero$RNE, NA_real_) ## used to be 0/0 = NaN

  ## a non-zero constant leaves NSE at rounding level rather than exactly zero,
  ## so RNE is meaningless there but must at least not be NaN
  one <- AdMitIS(N = 1e3, KERNEL = GelmanMeng, mit = fit$mit,
                 G = function(theta) rep(1, nrow(theta)))
  expect_equal(one$ghat, 1)
  expect_false(is.nan(one$RNE))
})

test_that("N is validated instead of being silently truncated", {
  mit <- list(p = 1, mu = as.matrix(0), Sigma = as.matrix(1), df = 1)
  expect_error(AdMitMH(N = 10.7, KERNEL = Gauss1d, mit = mit), "whole number")
  expect_error(AdMitIS(N = 10.7, KERNEL = Gauss1d, mit = mit), "whole number")
  expect_error(rMit(2.5, mit), "whole number")
  expect_error(AdMitMH(N = Inf, KERNEL = Gauss1d, mit = mit), "finite")
  expect_error(AdMitMH(N = c(10, 20), KERNEL = Gauss1d, mit = mit), "single number")
  expect_error(AdMit(KERNEL = GelmanMeng, mu0 = c(0, 0), control = list(Ns = 10000.5)),
               "whole number")
})

test_that("Sigma0 must be square and conform to mu0", {
  expect_error(AdMit(KERNEL = GelmanMeng, mu0 = c(0, 0), Sigma0 = matrix(1, 2, 3)),
               "must be a square matrix")
  expect_error(AdMit(KERNEL = GelmanMeng, mu0 = c(0, 0), Sigma0 = diag(3)),
               "as many rows as 'mu0'")
})

test_that("malformed mixtures are rejected with a useful message", {
  ok <- list(p = c(0.5, 0.5), mu = rbind(c(0, 0), c(1, 1)),
             Sigma = rbind(c(1, 0, 0, 1), c(1, 0, 0, 1)), df = 5)
  expect_silent(rMit(5, ok))

  bad <- ok; bad$mu <- ok$mu[1, , drop = FALSE]
  expect_error(dMit(c(0, 0), bad), "one row per mixture component")
  bad <- ok; bad$Sigma <- ok$Sigma[, 1:3]
  expect_error(dMit(c(0, 0), bad), "columns")
  bad <- ok; bad$df <- c(1, 2, 3)
  expect_error(dMit(c(0, 0), bad), "length 1 or 2")
  bad <- ok; bad$df <- -1
  expect_error(dMit(c(0, 0), bad), "finite and positive")
  bad <- ok; bad$p <- c(0.5, -0.5)
  expect_error(dMit(c(0, 0), bad), "finite and non-negative")
  bad <- ok; bad$mu <- as.vector(ok$mu)
  expect_error(dMit(c(0, 0), bad), "must be a matrix")
})