pkg/tests/testthat/test-jointDiag.R

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