tests/testthat/test-numerics.R

In AdMit: Adaptive Mixture of Student-t Distributions

## Regressions for the defects found in the 2.1.11 audit. All predate 2.1.10.

test_that("the fit is invariant to the scale of the log-kernel", {
  ## exp(2 * lnw) in fnlnf_C overflows past about +350 and underflows past
  ## about -350, so AdMit used to abort with "NA/NaN gradient evaluation" for
  ## any posterior with more than a handful of observations
  Shifted <- function(x, shift = 0, log = TRUE) {
    if (is.vector(x)) x <- matrix(x, nrow = 1)
    r <- -.5 * (x[,1]^2 * x[,2]^2 + x[,1]^2 + x[,2]^2 - 6 * x[,1] - 6 * x[,2]) + shift
    if (!log) r <- exp(r)
    as.vector(r)
  }
  set.seed(1234)
  base <- AdMit(Shifted, mu0 = c(0, 0.1), shift = 0,
                control = list(Ns = 5e3, Np = 5e2, Hmax = 3))
  for (s in c(-5000, -400, 400, 5000)) {
    set.seed(1234)
    out <- AdMit(Shifted, mu0 = c(0, 0.1), shift = s,
                 control = list(Ns = 5e3, Np = 5e2, Hmax = 3))
    ## the fit is close but not identical: the mode search in 'fn.optmu' uses
    ## optim()'s relative 'reltol', which is looser for a kernel offset by a
    ## few thousand. That is a property of the mode search, not of the
    ## probability optimizer tested below.
    expect_equal(out$CV, base$CV, tolerance = 0.02)
    expect_equal(out$mit$p, base$mit$p, tolerance = 0.02)
    expect_length(out$mit$p, length(base$mit$p))
  }
})

test_that("the probability objective is exactly invariant to a log-kernel shift", {
  set.seed(704)
  Np <- 200L; H <- 2L
  lnK <- matrix(rnorm(Np * H), Np, H)
  lnD <- matrix(rnorm(Np * H * H, -1, 0.5), Np, H * H)
  lnfgrad <- get("fn.lnfgrad", envir = asNamespace("AdMit"))

  base <- lnfgrad(lnK, lnD)
  for (s in c(-5000, -400, 400, 5000)) {
    obj <- lnfgrad(lnK + s, lnD)
    expect_equal(obj$f(0.3), base$f(0.3), tolerance = 1e-10)
    expect_equal(obj$g(0.3), base$g(0.3), tolerance = 1e-10)
    expect_true(is.finite(obj$f(0.3)))
  }
})

test_that("the objective and gradient survive extreme logits", {
  set.seed(701)
  Np <- 100L; H <- 2L
  lnK <- matrix(rnorm(Np * H), Np, H)
  lnD <- matrix(rnorm(Np * H * H, -1, 0.5), Np, H * H)
  obj <- get("fn.lnfgrad", envir = asNamespace("AdMit"))(lnK, lnD)

  ## exp(710) overflows: the softmax used to return c(NaN, 0)
  for (l in c(0, 100, 710, 1e4)) {
    expect_true(all(is.finite(obj$p(l))))
    expect_equal(sum(obj$p(l)), 1)
  }
  expect_true(is.finite(obj$f(0)))
  expect_true(all(is.finite(obj$g(0))))
})

test_that("a non-vectorized or malformed KERNEL is rejected", {
  mit <- list(p = 1, mu = matrix(0), Sigma = matrix(1), df = 1)
  ## a scalar result used to be recycled to length N and silently accepted
  expect_error(AdMitMH(500, function(theta, log = TRUE) 0, mit), "one value per point")
  ## a zero-length result left fnMH_C reading past the end of the vector
  expect_error(AdMitMH(500, function(theta, log = TRUE) numeric(0), mit), "one value per point")
  expect_error(AdMitIS(500, function(theta, log = TRUE) numeric(0), mit = mit), "one value per point")
  expect_error(AdMitMH(500, function(theta, log = TRUE) rep(NA_real_, nrow(as.matrix(theta))), mit),
               "NA or NaN")
  ## -Inf is a legitimate zero density and must still be allowed through
  expect_no_error(AdMitMH(200, function(theta, log = TRUE) {
    r <- -0.5 * as.matrix(theta)[,1]^2; r[1] <- -Inf; r }, mit))
})

test_that("a malformed G is rejected with an informative message", {
  set.seed(702)
  fit <- AdMit(KERNEL = GelmanMeng, mu0 = c(0, 0.1),
               control = list(Ns = 5e3, Np = 5e2, Hmax = 2))
  expect_error(AdMitIS(1000, GelmanMeng, mit = fit$mit, G = function(theta) 1),
               "one value per draw")
  expect_error(AdMitIS(1000, GelmanMeng, mit = fit$mit,
                       G = function(theta) rep(1, nrow(theta) - 1)), "one value per draw")
  expect_error(AdMitIS(1000, GelmanMeng, mit = fit$mit,
                       G = function(theta) rep(NA_real_, nrow(theta))), "NA or NaN")
  expect_error(AdMitIS(1000, GelmanMeng, mit = fit$mit, G = function(theta) "a"),
               "numeric vector or matrix")
})

test_that("mixture probabilities describe the same object in dMit and rMit", {
  ## dMit used the weights as given while rMit normalised them through sample()
  m <- list(p = c(2, 3), mu = matrix(c(0, 0), 2, 1),
            Sigma = matrix(c(1, 1), 2, 1), df = 5)
  expect_warning(dMit(0, m), "does not sum to one")
  ## integrate() calls dMit many times, each warning again
  v <- suppressWarnings(integrate(function(x) dMit(x, m, log = FALSE), -Inf, Inf)$value)
  expect_equal(v, 1, tolerance = 1e-6)
  ## a mixture that already sums to one is left strictly alone
  ok <- list(p = c(0.4, 0.6), mu = matrix(c(0, 0), 2, 1),
             Sigma = matrix(c(1, 1), 2, 1), df = 5)
  expect_no_warning(dMit(0, ok))
})

test_that("invalid scale matrices are rejected rather than giving Inf", {
  expect_error(dMit(0, list(p = 1, mu = matrix(0), Sigma = matrix(-1), df = 5)),
               "not positive definite")
  expect_error(rMit(3, list(p = 1, mu = matrix(0), Sigma = matrix(-1), df = 5)),
               "not positive definite")
  expect_error(dMit(c(0, 0), list(p = 1, mu = matrix(c(0, 0), 1, 2),
                                  Sigma = matrix(c(1, 1, 1, 1), 1, 4), df = 5)),
               "not positive definite")
  expect_error(dMit(c(0, 0), list(p = 1, mu = matrix(c(0, 0), 1, 2),
                                  Sigma = matrix(c(1, 2, 0, 1), 1, 4), df = 5)),
               "not symmetric")
  expect_error(dMit(0, list(p = 1, mu = matrix(NA_real_), Sigma = matrix(1), df = 5)),
               "finite and numeric")
})

test_that("the acceptance rate counts proposals, not draws", {
  ## the chain makes N-1 proposals; the rate used to be divided by N, so a
  ## perfect candidate could never report 1
  mit <- list(p = 1, mu = matrix(0), Sigma = matrix(1), df = 1)
  Ktarget <- function(x, log = TRUE) dMit(x, mit, log = log)
  for (n in c(2L, 10L, 100L)) {
    set.seed(1)
    expect_equal(AdMitMH(n, Ktarget, mit)$accept, 1)
  }
})

test_that("an exact candidate gives a zero coefficient of variation", {
  ## fn.CV used to stop with "'w' is constant"
  cv <- get("fn.CV", envir = asNamespace("AdMit"))
  expect_equal(cv(rep(2, 10)), 0)
  expect_error(cv(rep(0, 10)), "zero")

  mit <- list(p = 1, mu = matrix(0), Sigma = matrix(1), df = 1)
  Ktarget <- function(x, log = TRUE) dMit(x, mit, log = log)
  set.seed(703)
  out <- AdMit(Ktarget, mu0 = 0, Sigma0 = matrix(1),
               control = list(Ns = 1e3, Np = 5e2, Hmax = 4))
  expect_equal(out$CV[1], 0)
  expect_length(out$mit$p, 1) ## nothing to improve, so no second component
})

test_that("the ARCH kernel of the demo does not underflow", {
  theta <- matrix(c(1, 4, 0, .5), 1, 4)
  prior <- matrix(c(1, 0), 1, 2)
  y <- c(0, 100)
  got <- .C("fnKernelMixtureArch_C", theta = as.double(as.vector(t(theta))), N = 1L,
            y = as.double(y), n = 2L, prior = as.double(as.vector(t(prior))),
            d = vector("double", 1), PACKAGE = "AdMit", NAOK = TRUE)$d
  t1 <- log(0.5) + (-0.5 * (y[2] / sqrt(1))^2 - log(sqrt(1)))
  t2 <- log(0.5) + (-0.5 * (y[2] / sqrt(4))^2 - log(sqrt(4)))
  m <- max(t1, t2)
  expect_true(is.finite(got))
  expect_equal(got, m + log(exp(t1 - m) + exp(t2 - m)))
})

test_that("N is bounded to the integer range and AdMitIS needs two draws", {
  mit <- list(p = 1, mu = matrix(0), Sigma = matrix(1), df = 1)
  expect_error(AdMitMH(3e9, Gauss1d, mit), "must not exceed")
  expect_error(AdMitIS(1, Gauss1d, mit = mit), "should be larger than 1")
})

test_that("AdMit rejects non-finite mu0 and a zero ISscale", {
  expect_error(AdMit(GelmanMeng, mu0 = c(0, NA)), "finite and numeric")
  expect_error(AdMit(GelmanMeng, mu0 = c(0, Inf)), "finite and numeric")
  expect_error(AdMit(GelmanMeng, mu0 = c(0, 0.1),
                     control = list(Ns = 1e3, Np = 5e2, ISscale = c(1, 0))),
               "must be positive")
})