Nothing
senv <- function(X, Y, u, asy = TRUE, init = NULL){
X <- as.matrix(X)
Y <- as.matrix(Y)
a <- dim(Y)
n <- a[1]
r <- a[2]
p <- ncol(X)
if (a[1] != nrow(X))
stop("X and Y should have the same number of observations.")
if (u > r || u < 0)
stop("u must be an integer between 0 and r.")
if (sum(duplicated(cbind(X, Y), MARGIN = 2)) > 0)
stop("Some responses also appear in the predictors, or there maybe duplicated columns in X or Y.")
if (!is.null(init)) {
if (nrow(init) != r || ncol(init) != u) stop("The dimension of init is wrong.")
}
sigY <- stats::cov(Y) * (n - 1)/n
sigYX <- stats::cov(Y, X) * (n - 1)/n
sigX <- stats::cov(X) * (n - 1)/n
invsigX <- chol2inv(chol(sigX))
betaOLS <- sigYX %*% invsigX
U <- tcrossprod(betaOLS, sigYX)
M <- sigY - U
covMatrix <- NULL
asySE <- NULL
ratio <- NULL
if (u == 0) {
Gammahat <- NULL
Gamma0hat <- diag(r)
etahat <- NULL
Omegahat <- NULL
Omega0hat <- sigY
Lambdahat <- diag(r)
muhat <- colMeans(Y)
betahat <- matrix(0, r, p)
Sigmahat <- sigY
tmp.MU <- eigen(sigY)
objfun <- sum(log(tmp.MU$values))
loglik <- -n * r/2 * (log(2 * pi) + 1) - n/2 * objfun
if (asy == T)
ratio <- matrix(1, r, p)
}
else if (u >= (r - (r-1)/p)) {
Gammahat <- diag(r)
Gamma0hat <- NULL
etahat <- betaOLS
Omegahat <- M
Omega0hat <- NULL
Lambdahat <- diag(r)
muhat <- colMeans(Y) - betaOLS %*% colMeans(X)
betahat <- betaOLS
Sigmahat <- M
tmp.M <- eigen(M)
objfun <- sum(log(tmp.M$values))
loglik <- -n * r/2 * (log(2 * pi) + 1) - n/2 * objfun
if (asy == T) {
covMatrix <- kronecker(invsigX, M)
asySE <- matrix(sqrt(diag(covMatrix)), nrow = r)
ratio <- matrix(1, r, p)
}
}
else {
tmp <- senvMU(X, Y, u, initial = init)
Gammahat <- tmp$Gammahat
Gamma0hat <- tmp$Gamma0hat
Lambdahat <- tmp$Lambda
d1 <- diag(Lambdahat)
invLambda <- diag(1 / d1)
etahat <- crossprod(Gammahat, invLambda) %*% betaOLS
betahat <- Lambdahat %*% Gammahat %*% etahat
muhat <- colMeans(Y) - betahat %*% colMeans(X)
M1 <- invLambda %*% M %*% invLambda
Omegahat <- crossprod(Gammahat, M1) %*% Gammahat
M2 <- invLambda %*% sigY %*% invLambda
Omega0hat <- crossprod(Gamma0hat, M2) %*% Gamma0hat
E1 <- Lambdahat %*% Gammahat
Sigma1 <- E1 %*% tcrossprod(Omegahat, E1)
E2 <- Lambdahat %*% Gamma0hat
Sigmahat <- Sigma1 + E2 %*% tcrossprod(Omega0hat, E2)
tmp.MU <- eigen(sigY)
objfun <- tmp$objfun + sum(log(tmp.MU$values))
loglik <- -n * r/2 * (log(2 * pi) + 1) - n/2 * objfun
if (asy == T) {
covMatrix <- kronecker(invsigX, M)
asyFm <- matrix(sqrt(diag(covMatrix)), nrow = r)
SigmaL <- Gammahat %*% tcrossprod(Omegahat,
Gammahat) + Gamma0hat %*% tcrossprod(Omega0hat, Gamma0hat)
invsig <- chol2inv(chol(Sigmahat))
Cr <- contr(r)
Er <- expan(r)
Eu <- expan(u)
Eru <- expan(r - u)
g11 <- kronecker(diag(p), Gammahat)
g12 <- kronecker(t(etahat), diag(r))
g221 <- kronecker(Gammahat %*% Omegahat, diag(r))
g222 <- kronecker(Gammahat, Gamma0hat %*%
tcrossprod(Omega0hat, Gamma0hat))
g22 <- 2 * Cr %*% (g221 - g222)
g23 <- Cr %*% kronecker(Gammahat, Gammahat) %*% Eu
g24 <- Cr %*% kronecker(Gamma0hat, Gamma0hat) %*% Eru
G_o <- cbind(rbind(g11, matrix(0, nrow = nrow(g22), ncol = ncol(g11))),
rbind(g12, g22),
rbind(matrix(0, nrow = nrow(g11), ncol = ncol(g23)), g23),
rbind(matrix(0, nrow = nrow(g11), ncol = ncol(g24)), g24))
h_o1 <- kronecker(Gammahat %*% etahat, diag(r))
h_o2 <- 2 * tcrossprod(kronecker(SigmaL, diag(r)), Cr)
th_o <- cbind(h_o1, h_o2)
h_o <- t(th_o)
d11 <- kronecker(diag(p), Lambdahat)
d22 <- Cr %*% kronecker(Lambdahat, Lambdahat) %*% Er
DL <- cbind(rbind(d11, matrix(0, nrow = nrow(d22), ncol = ncol(d11))),
rbind(matrix(0, nrow = nrow(d11), ncol = ncol(d22)), d22))
diagvec <- function(d){
E <- NULL
for (i in 1 : (d - 1)){
e <- rep(0, d)
e[i + 1] <- 1
t <- kronecker(e, e)
E <- cbind(E, t)
}
E <- matrix(E, nrow = d^2)
}
L <- diagvec(r)
h1 <- DL %*% h_o %*% kronecker(diag(r), invLambda) %*% L
h2 <- DL %*% G_o
H <- cbind(h1, h2)
j11 <- kronecker(sigX, invsig)
j221 <- kronecker(invsig, invsig)
j22 <- 1/2 * crossprod(Er, j221) %*% Er
J <- cbind(rbind(j11, matrix(0, nrow = nrow(j22), ncol = ncol(j11))),
rbind(matrix(0, nrow = nrow(j11), ncol = ncol(j22)), j22))
M1 <- crossprod(H, J) %*% H
invM1 <- ginv(M1)
V1 <- tcrossprod(invM1, H)
V <- H %*% V1
covMatrix <- V[1:(p*r), 1:(p*r)]
asySE <- matrix(sqrt(diag(covMatrix)), nrow = r)
ratio <- asyFm/asySE
}
}
return(list(beta = betahat, Sigma = Sigmahat, Lambda = Lambdahat,
Gamma = Gammahat, Gamma0 = Gamma0hat, eta = etahat,
Omega = Omegahat, Omega0 = Omega0hat, mu = muhat,
loglik = loglik, covMatrix = covMatrix,
asySE = asySE, ratio = ratio, n = n))
}
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.