tests/testthat/test-MatrixNormal.R

In mniw: The Matrix-Normal Inverse-Wishart Distribution

#require(testthat)
## library(mniw)
source("mniw-testfunctions.R")
context("Matrix Normal Distribution")

tol <- 1e-6

test_that("Matrix Normal density is same in C++ as R", {
  calc.diff <- FALSE
  case.par <- expand.grid(p = c(1,2,4), q = c(1,2,3),
                          X = c("single", "multi"),
                          Lambda = c("none", "single", "multi"),
                          SigmaR = c("none", "single", "multi"),
                          SigmaC = c("none", "single", "multi"),
                          drop = c(TRUE, FALSE), stringsAsFactors = FALSE)
  ncases <- nrow(case.par)
  n <- 10
  if(calc.diff) {
    MaxDiff <- rep(NA, ncases)
  }
  for(ii in 1:ncases) {
    cp <- case.par[ii,]
    p <- cp$p
    q <- cp$q
    args <- list(X = list(p = p, q = q, rtype = cp$X, vtype = "matrix"),
                 Lambda = list(p = p, q = q, rtype = cp$Lambda, vtype = "matrix"),
                 SigmaR = list(p = p, rtype = cp$SigmaR, vtype = "matrix"),
                 SigmaC = list(q = q, rtype = cp$SigmaC, vtype = "matrix"))
    args <- get_args(n = n, args = args, drop = cp$drop)
    # R test
    llR <- rep(NA, n)
    for(jj in 1:n) {
      llR[jj] <- dMNormR(X = args$R$X[[jj]],
                         Lambda = args$R$Lambda[[jj]],
                         SigmaU = args$R$SigmaR[[jj]],
                         SigmaV = args$R$SigmaC[[jj]],
                         log = TRUE)
    }
    # C++ test
    llcpp <- do.call(dMNorm, args = c(args$cpp, list(log = TRUE)))
    # if all inputs to c++ are single it returns only one value
    if(all_single(cp)) {
      llcpp <- rep(llcpp, n)
    }
    mx <- abs(llR-llcpp)
    mx <- min(max(mx), max(mx/abs(llR)))
    if(calc.diff) {
      MaxDiff[ii] <- mx
    } else {
      ## expect_equal(mx, 0, tolerance = tol)
      expect_Rcpp_equal("dMNorm", ii, mx, tolerance = tol)
    }
  }
})

test_that("Matrix Normal simulation is same in C++ as R", {
  calc.diff <- FALSE
  case.par <- expand.grid(p = c(1,2,4), q = c(1,2,3),
                          Lambda = c("none", "single", "multi"),
                          SigmaR = c("none", "single", "multi"),
                          SigmaC = c("none", "single", "multi"),
                          drop = c(TRUE, FALSE), stringsAsFactors = FALSE)
  # remove cases where dimensions can't be identified
  case.par <- case.par[!with(case.par, {
    Lambda == "none" & ((SigmaR == "none") | (SigmaC == "none"))}),]
  ncases <- nrow(case.par)
  rownames(case.par) <- 1:ncases
  n <- 10
  TestSeed <- sample(1e6, ncases)
  if(calc.diff) {
    MaxDiff <- rep(NA, ncases)
  }
  for(ii in 1:nrow(case.par)) {
    set.seed(TestSeed[ii]) # seed
    cp <- case.par[ii,]
    p <- cp$p
    q <- cp$q
    args <- list(Lambda = list(p = p, q = q, rtype = cp$Lambda, vtype = "matrix"),
                 SigmaR = list(p = p, rtype = cp$SigmaR, vtype = "matrix"),
                 SigmaC = list(q = q, rtype = cp$SigmaC, vtype = "matrix"))
    args <- get_args(n = n, args = args, drop = cp$drop)
    # R test
    set.seed(TestSeed[ii]) # seed
    XR <- array(NA, dim = c(p,q,n))
    for(jj in 1:n) {
      XR[,,jj] <- rMNormR(Lambda = args$R$Lambda[[jj]],
                          SigmaU = args$R$SigmaR[[jj]],
                          SigmaV = args$R$SigmaC[[jj]])
    }
    # C++ test
    set.seed(TestSeed[ii]) # seed
    Xcpp <- do.call(rMNorm, args = c(args$cpp, list(n = n)))
    mx <- abs(range(XR - Xcpp))
    mx <- min(max(mx), max(mx/abs(XR)))
    if(calc.diff) {
      MaxDiff[ii] <- mx
    } else {
      ## expect_equal(mx, 0, tolerance = tol)
      expect_Rcpp_equal("rMNorm", ii, mx, tolerance = tol)
    }
  }
})