R/deep.sem.alg.2.R

Defines functions deep.sem.alg.2

deep.sem.alg.2 <- function(y, numobs, p, r, k, H.list, psi.list, psi.list.inv,
                           mu.list, w.list, it, eps) {

likelihood <- NULL
hh <- 0
ratio <- Inf #1000
layers <- length(k)
#################################
#### compute the likelihood #####
#################################
out <- compute.lik(y, numobs, k, mu.list, H.list, psi.list, w.list)
py <- out$py
ps.y <- out$ps.y
ps.y.list <- out$ps.y.list
k.comb <- out$k.comb
s <- out$s
tot.k <- prod(k)
temp <- sum(log(py))
likelihood <- c(likelihood, temp)

z.list <- NULL
#####################################################
while ((hh < it) & (ratio > eps )) {
  hh <- hh + 1
  ###############################################################
  ###################  first layer ##############################
  ###############################################################
  l <- 1
  yy <- y
  z2 <- z2.one <- array(0, c(numobs, r[l + 1], k[l], k[l + 1]))
  zz2 <- array(0, c(numobs, r[l + 1], r[l + 1], k[l], k[l + 1]))
  z <- z.one <- array(0, c(numobs, r[l + 1], k[l]))
  zz <- array(0, c(numobs, r[l + 1], r[l + 1], k[l]))

  for (p1 in 1 : k[l]) {
    for (p2 in 1 : k[l + 1]) {

      # A is the inverse of \xi_{sl}
      # = (\Sigma_{sl}^(-1) + \Lambda_{sl}^T \Psi_{sl}^(-1) \Lambda_{sl} )
      A <- ginv(H.list[[l + 1]][p2,, ] %*% t(H.list[[l + 1]][p2,, ]) +
            psi.list[[l + 1]][p2,, ]) + t(H.list[[l]][p1,, ]) %*%
            (psi.list.inv[[l]][p1,, ]) %*% (H.list[[l]][p1,, ])

      # b contains part of p_{sl} (z^{l-1})
      # =  Sigma_{sl}^(-1) mu_{sl}
      b <- ginv(H.list[[l + 1]][p2,, ] %*% t(H.list[[l + 1]][p2,, ]) +
            psi.list[[l + 1]][p2,, ]) %*%
            matrix(mu.list[[l + 1]][, p2], r[l + 1], numobs) +
            t(H.list[[l]][p1,, ]) %*% (psi.list.inv[[l]][p1,, ]) %*%
            (t(yy) - matrix(mu.list[[l]][, p1], r[l], numobs))

      chsi <- ginv(A)
      if (!isSymmetric(chsi)) {
        chsi <- makeSymm(chsi)
      }
      roy <- chsi %*% b
      roy.quadro <- array(apply(roy, 2, function(x) x %*% t(x)),
                         c(r[l + 1], r[l + 1], numobs))
      zz2[,,, p1, p2] <- aperm(array(chsi, c(r[l + 1], r[l + 1], numobs)) +
                         roy.quadro, c(3, 1, 2))

      if (!is.positive.definite(chsi)) {
        chsi <- make.positive.definite(chsi)
      }

      z2.one[,, p1, p2] <- rmvnorm(numobs, rep(0, r[l + 1]), chsi) + t(roy)
      z2[,, p1, p2] <- t(roy)
    }
  }

for (i in 1 : k[l + 1]) {
  prob <- ps.y.list[[l + 1]][, `i`, drop = FALSE]
  z <- z + array(z2[,,, i, drop = FALSE] *
           array(rowSums(prob), c(numobs, r[l + 1], k[l], 1)),
           c(numobs, r[l + 1], k[l]))
  z.one <- z.one + array(z2.one[,,, i, drop = FALSE] *
                   array(rowSums(prob), c(numobs, r[l + 1], k[l], 1)),
                   c(numobs, r[l+1], k[l]))
  zz <- zz + array(zz2[,,,, i, drop = FALSE] *
             array(rowSums(prob), c(numobs, r[l + 1], r[l + 1], k[l], 1)),
            c(numobs, r[l + 1], r[l + 1], k[l]))
}

z.list[[l]] <- aperm(z.one, c(3, 1, 2))

out <- compute.est(k[l], r[l], r[l+1], ps.y.list[[l]], yy,
        aperm(z, c(3, 1, 2)), aperm(zz, c(4, 2, 3, 1)), mu.list[[l]])

H.list[[l]] <- out$H
psi.list[[l]] <- out$psi
psi.list.inv[[l]] <- out$psi.inv
mu.list[[l]] <- out$mu
w.list[[l]] <- out$w

###############################################################
################### second layer ##############################
###############################################################
l <- 2
yy <- matrix(0, numobs, r[l])
zz <- z.list[[l - 1]]
for (i in 1 : k[l - 1]) {
  yy <- yy + matrix(zz[i, ,, drop = FALSE] *
           array(rowSums(ps.y.list[[l - 1]][, i, drop = FALSE]),
           c(1, numobs, r[l])), numobs, r[l])
}

z <- z.one <- array(0, c(numobs, r[l + 1], k[l]))
zz <- array(0, c(numobs, r[l + 1], r[l + 1], k[l]))

for (p1 in 1 : k[l]) {
  A <- diag(r[l + 1]) + t(H.list[[l]][p1,, ]) %*%
       (psi.list.inv[[l]][p1,, ]) %*% (H.list[[l]][p1,, ])
  b <- t(H.list[[l]][p1,, ]) %*% (psi.list.inv[[l]][p1,, ]) %*%
       (t(yy) - matrix(mu.list[[l]][, p1], r[l], numobs))

  chsi <- ginv(A)

  if (!isSymmetric(chsi)) {
    chsi <- makeSymm(chsi)
  }

  roy <- chsi %*% b
  roy.quadro <- array(apply(roy, 2, function(x) x %*% t(x)),
                  c(r[l + 1], r[l + 1], numobs))
  zz[,,,p1] <- aperm(array(chsi, c(r[l + 1], r[l + 1], numobs)) +
                 roy.quadro,c(3, 1, 2))

  if (!is.positive.definite(chsi)) {
      chsi <- make.positive.definite(chsi)
  }

  z.one[,, p1] <- t(roy) + rmvnorm(numobs, rep(0, r[l + 1]), chsi)
  z[,, p1] <- t(roy)
}

z.list[[l]] <- aperm(z.one, c(3, 1, 2))
out <- compute.est(k[l], r[l], r[l + 1], ps.y.list[[l]], yy,
        aperm(z, c(3, 1, 2)), aperm(zz, c(4, 2, 3, 1)), mu.list[[l]])

H.list[[l]] <- out$H
psi.list[[l]] <- out$psi
psi.list.inv[[l]] <- out$psi.inv
mu.list[[l]] <- out$mu
w.list[[l]] <- out$w

out <- compute.lik(y, numobs, k, mu.list, H.list, psi.list, w.list)
py <- out$py
ps.y <- out$ps.y
ps.y.list <- out$ps.y.list
k.comb <- out$k.comb
s <- out$s

lik <- sum(log(py))
likelihood <- c(likelihood, lik)

if (hh < 5) {
  ratio <- 2 * eps
}
if (hh > 5) {
  ratio <- (ma(likelihood, 5) [hh + 1] - ma(likelihood, 5) [hh]) /
            abs(ma(likelihood,5)[hh])
}

}

h <- 0
for (j in 1 : layers) {
  h <- h + (k[j] - 1) + (r[j] * r[j + 1]) * k[j] + r[j] * k[j] + k[j] * r[j]
}
if (layers > 1) {
  for (j in 2:layers) {
    h <- h - (r[j] * k[j] * (r[j] - 1) / 2)
  }
}

bic <- -2 * lik + h * log(numobs)
aic <- -2 * lik + 2 * h
EN <- entr(ps.y.list[[1]])
clc <- -2 * lik + 2 * EN
icl.bic <- -2 * lik + 2 * EN + h * log(numobs)

out <- list(H = H.list, w = w.list, mu = mu.list, psi = psi.list,
            likelihood = likelihood, bic = bic, aic = aic, clc = clc,
            s = s, icl_bic = icl.bic, h = h, ps.y = ps.y)
return(out)
}

Try the deepgmm package in your browser

Any scripts or data that you put into this service are public.

deepgmm documentation built on Nov. 21, 2022, 1:05 a.m.