Nothing
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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.