## 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") })