pkg/tests/testthat/test-jointDiag.R

   1 #auxiliary to test diagonality
   2 .computeMuCheckDiag = function(X, Y, K, jd_method, β_ref)
   3 {
   4   d <- ncol(X)
   5   #TODO: redundant code, same as computeMu() main method. Comments are stripped
   6   M3 <- .Moments_M3(X,Y)
   7   M2_t <- array(dim=c(d,d,d))
   8   for (i in 1:d)
   9   {
  10     ρ <- c(rep(0,i-1),1,rep(0,d-i))
  11     M2_t[,,i] <- .T_I_I_w(M3,ρ)
  12   }
  13   jd <- jointDiag::ajd(M2_t, method=jd_method)
  14   V <- if (jd_method=="uwedge") jd$B else solve(jd$A)
  15   M2_t <- array(dim=c(d,d,K))
  16   for (i in 1:K)
  17     M2_t[,,i] <- .T_I_I_w(M3,V[,i])
  18   #END of computeMu() code
  19
  20   max_error <- 1.0 #TODO: tune ?
  21   invβ <- MASS::ginv(β_ref)
  22   for (i in 1:K)
  23   {
  24     shouldBeDiag <- invβ %*% M2_t[,,i] %*% t(invβ)
  25     expect_lt(
  26       mean( abs(shouldBeDiag[upper.tri(shouldBeDiag) | lower.tri(shouldBeDiag)]) ),
  27       max_error)
  28   }
  29 }
  30
  31 test_that("'jedi' and 'uwedge' joint-diagonalization methods return a correct matrix",
  32 {
  33   n <- 10000
  34   d_K <- list( c(2,2), c(5,3), c(20,13) )
  35
  36   for (dk_index in 1:length(d_K))
  37   {
  38     d <- d_K[[dk_index]][1]
  39     K <- d_K[[dk_index]][2]
  40     #NOTE: sometimes large errors if pr is not balanced enough (e.g. random);
  41     #      same note for β. However we could be more random than that...
  42     β_ref <- rbind(diag(K),matrix(0,nrow=d-K,ncol=K))
  43     io <- generateSampleIO(n, p=rep(1/K,K-1), β=β_ref, rep(0,K), link="logit")
  44     .computeMuCheckDiag(io$X, io$Y, K, jd_method="uwedge", β_ref)
  45     #TODO: some issues with jedi method (singular system)
  46     #.computeMuCheckDiag(io$X, io$Y, K, jd_method="jedi", β_ref)
  47   }
  48 })