mit1d <- list(p = c(0.5, 0.5),
              mu = matrix(c(-2.0, 0.5), 2, 1, byrow = TRUE),
              Sigma = matrix(0.1, 2, 1),
              df = 10)

mit2d <- 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)

test_that("rMit returns an Nxk matrix for a multivariate mixture", {
  set.seed(201)
  expect_identical(dim(rMit(10, mit2d)), c(10L, 2L))
  ## regression: a single draw used to collapse to a length-k vector
  expect_identical(dim(rMit(1, mit2d)), c(1L, 2L))
  expect_identical(dim(rMit(2, mit2d)), c(2L, 2L))
})

test_that("rMit returns a vector of length N for a univariate mixture", {
  set.seed(202)
  expect_null(dim(rMit(10, mit1d)))
  expect_length(rMit(10, mit1d), 10)
  expect_length(rMit(1, mit1d), 1)
})

test_that("rMit draws match the mixture moments", {
  set.seed(203)
  theta <- rMit(2e4, mit2d)
  expected_mean <- colSums(mit2d$p * mit2d$mu)
  expect_equal(colMeans(theta), expected_mean, tolerance = 0.05,
               ignore_attr = TRUE)
})

test_that("rMit honours component-specific degrees of freedom", {
  set.seed(204)
  ## df supplied as a vector of length H, as allowed since version 1-01.04
  mit <- mit2d
  mit$df <- c(3, 30)
  expect_identical(dim(rMit(50, mit)), c(50L, 2L))
})

test_that("dMit matches the mixture density by construction", {
  theta <- rbind(c(0, 0), c(2, -1), c(1.5, 0.5))
  ref <- rowSums(vapply(seq_along(mit2d$p), function(h)
    mit2d$p[h] * mvtnorm::dmvt(theta, mit2d$mu[h, ],
                               matrix(mit2d$Sigma[h, ], 2, 2),
                               mit2d$df, log = FALSE), numeric(nrow(theta))))

  expect_equal(dMit(theta, mit2d, log = FALSE), ref)
  expect_equal(dMit(theta, mit2d, log = TRUE), log(ref))
})

test_that("dMit accepts a single point supplied as a vector", {
  expect_equal(dMit(c(2, -1), mit2d), dMit(matrix(c(2, -1), 1), mit2d))
  expect_length(dMit(c(2, -1), mit2d), 1)
})

test_that("dMit integrates to one over the univariate mixture", {
  f <- function(x) dMit(x, mit1d, log = FALSE)
  expect_equal(integrate(f, -Inf, Inf)$value, 1, tolerance = 1e-4)
})

test_that("rMit and dMit warn and fall back when mit is empty", {
  set.seed(205)
  expect_warning(rMit(5, list()), "not well defined")
  expect_warning(dMit(0, list()), "not well defined")
})