# R/genv.R In Renvlp: Computing Envelope Estimators

#### Documented in genv

```genv <- function(X, Y, Z, u, asy = TRUE, fit = TRUE, init = NULL){

groupind <- unique(Z)
ZZ <- as.factor(Z)
Z <- match(ZZ, levels(ZZ))
ZZ <- as.factor(Z)
X <- as.matrix(X)
Y <- as.matrix(Y)
a <- dim(Y)
n <- a[1]
r <- a[2]
c <- ncol(X)
p <- nlevels(ZZ)
if (!is.null(init)) {
if (nrow(init) != r || ncol(init) != u) stop("The dimension of init is wrong.")
}

if(u < 0 | u > r){
print("u should be an interger between 0 and r")
skip <- 1
} else {
skip <- 0
}
if (n <= c) {
print("This works when p < n")
skip <- 1
} else {
skip <-0
}

if(skip == 0){

ncumx <- c()
for (i in 1 : p) {
ncumx[i] <- length(which(ZZ == as.numeric(levels(ZZ)[i])))
}
ncum <- cumsum(ncumx)
ng <- diff(c(0, ncum))
sortz <- sort(Z, index.return = T)
Zs <- sortz\$x
ind <- sortz\$ix
Ys <- Y[ind, ]
Xs <- X[ind, ]
mY <- apply(Y, 2, mean)
zz <- list(length = p)
for(i in 1:p){
if(i > 1) {
zz[[i]] <- stats::cov(Xs[((ncum[i - 1] + 1):ncum[i]), ])
}else{
zz[[i]] <- stats::cov(Xs[(1:ncum[i]), ])
}
}
sigres = sigYcg = sigresc <- list(length = p)
mYg <- matrix(rep(0, r*p), ncol = p)
mXg <- matrix(rep(0, c*p), ncol = p)
for (i in 1 : p) {
if(i > 1) {
yy <- Ys[(ncum[i-1] + 1):ncum[i], ]
xx <- Xs[(ncum[i-1] + 1):ncum[i], ]
sigres[[i]] <- stats::cov(yy) * (ng[i] - 1)/ng[i]
mYg[ , i] <- apply(yy, 2, mean)
mXg[ , i] <- apply(xx, 2, mean)
Ycg <- scale(yy, center=T, scale=F)
Xcg <- scale(xx, center=T, scale=F)
QZC <- diag(ng[i]) - Xcg %*% solve(t(Xcg) %*% Xcg) %*% t(Xcg)
sigYcg[[i]] <- t(Ycg) %*% Ycg
sigresc[[i]] <- t(Ycg) %*% QZC %*% Ycg/ng[i]
} else {
yy <- Ys[1:ncum[i], ]
xx <- Xs[1:ncum[i], ]
sigres[[i]]<- stats::cov(yy) * (ng[i] - 1)/ng[i]
mYg[ , i] <- apply(yy, 2, mean)
mXg[ , i] <- apply(xx, 2, mean)
Ycg <- scale(yy, center=T, scale=F)
Xcg <- scale(xx, center=T, scale=F)
QZC <- diag(ng[i]) - Xcg %*% solve(t(Xcg) %*% Xcg) %*% t(Xcg)
sigYcg[[i]] <- t(Ycg) %*% Ycg
sigresc[[i]] <- t(Ycg) %*% QZC %*% Ycg/ng[i]
}
}
sigY <- stats::cov(Y) * (n - 1)/n
eigtemY <- eigen(sigY)\$values
logDetSigY <- log(prod(eigtemY[eigtemY > 0]))
invsigY <- solve(sigY)

sigYc <- matrix(rep(0, r*r), ncol = r)
for (i in 1 : p){
sigYc <- sigYc + sigYcg[[i]]/n
}
eigtemYc <- eigen(sigYc)\$values
logDetSigYc <- log(prod(eigtemYc[eigtemYc > 0]))
invsigYc <- solve(sigYc)

U = M <- list(length = p)
for (i in 1 : p){
M[[i]] <- sigresc[[i]]
U[[i]] <-  sigYc - M[[i]]
}
MU <- sigYc
tmp <- genvMU(M, U, MU, u, n, ng, p, initial = init)

Gammahat <- tmp\$Gammahat
Gamma0hat <- tmp\$Gamma0hat
covMatrix <- NULL
asySE <- NULL
ratio <- NULL
Yfit <- NULL

if(u == 0){

Sigmahat <- sigYc
etahat <- NULL
Omegahat <- NULL
Omega0hat <- sigYc
betahat <- list(length=p)
for (i in 1:p){
betahat[[i]] <- matrix(rep(0, r*c), ncol=c)
}
muhat <- mYg
bb <- rep(0, p)
for (i in 1:p){
if(i > 1){
yy <- Ys[(ncum[i - 1]+1):ncum[i],]
Ycg <- scale(yy, center = T, scale = F)
bb[i] <- sum(diag(Ycg %*% solve(sigYc) %*% t(Ycg)))
} else {
yy <- Ys[1 : ncum[i], ]
Ycg <- scale(yy, center = T, scale = F)
bb[i] <- sum(diag(Ycg %*% solve(sigYc) %*% t(Ycg)))
}
}
b <- sum(bb)
loglik <- - n * r * log(2 * pi)/2 - n * logDetSigYc/2 - b/2
if (asy == T) {
ratio <- list(length = p)
for(i in 1:p) {
ratio[[i]] <- matrix(1, r, c)
}
}
if (fit == T) {
Yfit <- matrix(rep(0, n * r), ncol = r)
}

} else if (u == r) {

Sigmahat <- sigresc
aa <- rep(0, p)
etahat <- list(length = p)
for (i in 1 : p){
if(i > 1){
yy <- Ys[(ncum[i - 1]+1):ncum[i], ]
xx <- Xs[(ncum[i - 1]+1):ncum[i], ]
Ycg <- scale(yy, center = T, scale = F)
Xcg <- scale(xx, center = T, scale = F)
etahat[[i]] <- t(Ycg) %*% Xcg %*% solve(t(Xcg) %*% Xcg)
QZC <- diag(ng[i]) - Xcg %*% solve(t(Xcg) %*% Xcg) %*% t(Xcg)
aa[i] <- sum(diag(QZC %*% Ycg %*%
solve(sigresc[[i]]) %*% t(Ycg) %*% QZC))
} else {
yy <- Ys[1:ncum[i], ]
xx <- Xs[1:ncum[i], ]
Ycg <- scale(yy, center = T, scale = F)
Xcg <- scale(xx, center = T, scale = F)
etahat[[i]] <- t(Ycg) %*% Xcg %*% solve(t(Xcg) %*% Xcg)
QZC <- diag(ng[i]) - Xcg %*% solve(t(Xcg) %*% Xcg) %*% t(Xcg)
aa[i] <- sum(diag(QZC %*% Ycg %*% solve(sigresc[[i]])
%*% t(Ycg) %*% QZC))
}
}
a <- sum(aa)
Omegahat <- sigresc
Omega0hat <- NULL
betahat <- etahat
muhat <- mYg
for (i in 1:p){
muhat[ , i] <- muhat[ , i] - betahat[[i]] %*% mXg[ , i]
}
loglik <- - n * r * log(2 * pi)/2 - a/2
for (i in 1 : p){
eig <- eigen(sigresc[[i]])
loglik <- loglik - ng[i]*sum(log(eig\$values))/2
}
if (asy == T) {
asybeta = sebeta = asySEbeta = ratio <- list(length = p)
for(i in 1:p){
asybeta[[i]] <- n * kronecker(solve(zz[[i]]), sigresc[[i]]) / ng[i]
sebeta[[i]] <- diag(asybeta[[i]])
asySEbeta[[i]] <- matrix(sqrt(sebeta[[i]]), ncol = c)
}
covMatrix <- asybeta
asySE <- asySEbeta
for(i in 1:p) {
ratio[[i]] <- matrix(1, r, c)
}
}
if (fit == T) {
Yfit <- matrix(rep(0, n * r), ncol = r)
for (i in 1 : p) {
if(i > 1) {
Yfit[ind[(ncum[i - 1] + 1):ncum[i]], ] <- rep(1, ng[i]) %*% t(muhat[ ,i]) + X[ind[(ncum[i - 1] + 1):ncum[i]],] %*% t(betahat[[i]])
} else {
Yfit[ind[1:ncum[1]], ] <- rep(1, ng[1]) %*% t(muhat[ ,1]) + X[ind[1:ncum[1]], ] %*% t(betahat[[1]])
}
}
}

} else {
Omega0hat <- t(Gamma0hat) %*% sigYc %*% Gamma0hat
Sigmahat = Omegahat <- list(length = p)
for(i in 1 : p){
Omegahat[[i]] <- t(Gammahat) %*% sigresc[[i]] %*% Gammahat
Sigmahat[[i]] <- Gammahat %*% Omegahat[[i]] %*% t(Gammahat) + Gamma0hat %*% Omega0hat %*% t(Gamma0hat)
}
aa = bb <- rep(0, p)
etahat = betahat <- list(length = p)
for (i in 1 : p){
if(i > 1){
yy <- Ys[(ncum[i - 1] + 1):ncum[i], ]
xx <- Xs[(ncum[i - 1] + 1):ncum[i], ]
Ycg <- scale(yy, center = T, scale = F)
Xcg <- scale(xx, center = T, scale = F)
etahat[[i]] <- t(Gammahat) %*% t(Ycg) %*% Xcg %*% solve(t(Xcg) %*% Xcg)
betahat[[i]] <- Gammahat %*% etahat[[i]]
QZC <- diag(ng[i]) - Xcg %*% solve(t(Xcg) %*% Xcg) %*% t(Xcg)
aa[i] <- sum(diag(QZC %*% Ycg %*% Gammahat %*% solve(t(Gammahat) %*% sigresc[[i]] %*%
Gammahat) %*% t(Gammahat) %*% t(Ycg) %*% QZC))
bb[i] <- sum(diag(Ycg %*% Gamma0hat %*% solve(t(Gamma0hat) %*%
sigYc %*% Gamma0hat) %*% t(Gamma0hat) %*% t(Ycg)))
}else{
yy <- Ys[1 : ncum[i], ]
xx <- Xs[1 : ncum[i], ]
Ycg <- scale(yy, center = T, scale = F)
Xcg <- scale(xx, center = T, scale = F)
etahat[[i]] <- t(Gammahat) %*% t(Ycg) %*% Xcg %*% solve(t(Xcg) %*% Xcg)
betahat[[i]] <- Gammahat %*% etahat[[i]]
QZC <- diag(ng[i]) - Xcg %*% solve(t(Xcg) %*% Xcg) %*% t(Xcg)
aa[i] <- sum(diag(QZC %*% Ycg %*% Gammahat %*%
solve(t(Gammahat) %*% sigresc[[i]] %*% Gammahat) %*% t(Gammahat) %*% t(Ycg) %*% QZC))
bb[i] <- sum(diag(Ycg %*% Gamma0hat %*% solve(t(Gamma0hat) %*%
sigYc %*% Gamma0hat) %*% t(Gamma0hat) %*% t(Ycg)))
}
}
muhat <- mYg
for (i in 1 : p){
muhat[ , i]<- muhat[ , i] - betahat[[i]] %*% mXg[ , i]
}
a <- sum(aa)
b <- sum(bb)
eig2 <- eigen(t(Gamma0hat) %*% sigYc %*% Gamma0hat)
r1 <- sum(log(eig2\$values))
loglik <- -n * r * log(2 * pi)/2 - a/2 -b/2 - n * r1/2
for (i in 1:p){
eig <- eigen(t(Gammahat) %*% sigresc[[i]] %*% Gammahat)
loglik <- loglik - ng[i] * sum(log(eig\$values))/2
}
if (asy == T) {
asybeta = sebeta = asyFmbeta <- list(length = p)
for(i in 1:p){
asybeta[[i]] <- n * kronecker(solve(zz[[i]]), sigresc[[i]]) / ng[i]
sebeta[[i]] <- diag(asybeta[[i]])
asyFmbeta[[i]] <- matrix(sqrt(sebeta[[i]]), ncol = c)
}
covMatrix <- asybeta
asyFm <- asyFmbeta

asybeta <- list(length = p)
sebeta = asySEbeta = ratio <- list(length = p)
bb <- matrix(rep(0, u * (r - u) * u * (r - u)), ncol = u * (r - u))
for(i in 1 : p) {
b <- kronecker(etahat[[i]] %*% zz[[i]] %*% t(etahat[[i]]),
solve(Omega0hat)) + kronecker(Omegahat[[i]],
solve(Omega0hat)) + kronecker(solve(Omegahat[[i]]),
Omega0hat) - 2 * kronecker(diag(u), diag(r - u))
bb <- bb + ng[i] * b/n
}
for(i in 1 : p) {
aux <- kronecker(t(etahat[[i]]),
Gamma0hat) %*% solve(bb) %*% kronecker(etahat[[i]],
t(Gamma0hat))
asybeta[[i]] <- n * kronecker(solve(zz[[i]]), Gammahat %*%
Omegahat[[i]] %*% t(Gammahat)) / ng[i] + aux
sebeta[[i]] <- diag(asybeta[[i]])
asySEbeta[[i]] <- matrix(sqrt(sebeta[[i]]), ncol = c)
}
covMatrix <- asybeta
asySE <- asySEbeta
for(i in 1:p){
ratio[[i]] <- asyFm[[i]] / asySE[[i]]
}
}
if (fit == T) {
Yfit <- matrix(rep(0, n*r), ncol = r)
for (i in 1 : p) {
if(i > 1) {
Yfit[ind[(ncum[i - 1] + 1):ncum[i]], ] <- rep(1, ng[i]) %*% t(muhat[ ,i]) + X[ind[(ncum[i - 1] + 1):ncum[i]], ] %*% t(betahat[[i]])
} else {
Yfit[ind[1:ncum[1]], ] <- rep(1, ng[1]) %*% t(muhat[ ,1]) + X[ind[1:ncum[1]], ] %*% t(betahat[[1]])
}
}
}

}

return(list(Gamma = Gammahat, Gamma0 = Gamma0hat,
beta = betahat, Sigma = Sigmahat,
eta = etahat, Omega = Omegahat, Omega0 = Omega0hat,
mu = muhat, loglik = loglik, n = n, ng = ng, Yfit = Yfit,
covMatrix = covMatrix, asySE = asySE, ratio = ratio, groupInd = groupind))
}
}
```

## Try the Renvlp package in your browser

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

Renvlp documentation built on Sept. 11, 2021, 9:07 a.m.