Nothing
stenv <- function (X, Y, q, u, asy = TRUE, Pinit = NULL, Ginit = 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 (q > p || q < 0)
stop("q must be an integer between 0 and p.")
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(Pinit)) {
if (nrow(Pinit) != r || ncol(Pinit) != u) stop("The dimension of Pinit is wrong.")
}
if (!is.null(Ginit)) {
if (nrow(Ginit) != p || ncol(Ginit) != q) stop("The dimension of Ginit 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
covMatrix <- NULL
asySE <- NULL
ratio <- NULL
if (u == 0 || q ==0) {
Gammahat <- NULL
Gamma0hat <- diag(r)
Phihat <- NULL
Phi0hat <- diag(p)
etahat <- NULL
Omegahat <- NULL
Omega0hat <- sigY
Deltahat <- NULL
Delta0hat <- sigX
muhat <- colMeans(Y)
betahat <- matrix(0, p, r)
SigmaYcXhat <- sigY
SigmaXhat <- sigX
tmp1 <- eigen(sigX)
tmp2 <- eigen(sigY)
objfun <- sum(log(tmp1$values)) + sum(log(tmp2$values))
loglik <- - n/2 * objfun
if (asy == T)
ratio <- matrix(1, p, r)
}
else if (u == r && q == p) {
M <- sigY - tcrossprod(betaOLS, sigYX)
Gammahat <- diag(r)
Gamma0hat <- NULL
Phihat <- diag(p)
Phi0hat <- NULL
etahat <- betaOLS
Omegahat <- M
Omega0hat <- NULL
Deltahat <- sigX
Delta0hat <- NULL
muhat <- colMeans(Y) - betaOLS %*% colMeans(X)
betahat <- t(betaOLS)
SigmaYcXhat <- M
SigmaXhat <- sigX
tmp.e1 <- eigen(Omegahat)
tmp.e3 <- eigen(Deltahat)
objfun <- sum(log(tmp.e1$values)) + sum(log(tmp.e3$values))
loglik <- -n * (p + r) / 2 - n/2 * objfun
if (asy == T) {
covMatrix <- kronecker(M, invsigX)
asySE <- matrix(sqrt(diag(covMatrix)), nrow = p)
ratio <- matrix(1, p, r)
}
}
else if (u == r) {
M1 <- sigY - tcrossprod(betaOLS, sigYX)
tmp.xenv <- xenv(X, Y, q, asy = TRUE, init = Ginit)
Gammahat <- diag(r)
Gamma0hat <- NULL
Phihat <- tmp.xenv$Gamma
Phi0hat <- tmp.xenv$Gamma0
E1 <- crossprod(Phihat, sigX) %*% Phihat
invE1 <- chol2inv(chol(E1))
etahat <- invE1 %*% t(Phihat) %*% t(sigYX)
Omegahat <- M1
Omega0hat <- NULL
Deltahat <- tmp.xenv$Omega
Delta0hat <- tmp.xenv$Omega0
muhat <- tmp.xenv$mu
betahat <- Phihat %*% etahat
SigmaYcXhat <- tmp.xenv$SigmaYcX
SigmaXhat <- tmp.xenv$SigmaX
invome <- chol2inv(chol(Omegahat))
invdel <- chol2inv(chol(Deltahat))
invdel0 <- chol2inv(chol(Delta0hat))
tmp.e1 <- eigen(Omegahat)
tmp.e2 <- eigen(Deltahat)
tmp.e3 <- eigen(Delta0hat)
t3 <- sigY %*% invome
t5 <- crossprod(betahat, sigX) %*% betahat %*% invome
t6 <- t(betahat) %*% t(sigYX) %*% invome
obj1 <- sum(log(tmp.e1$values)) + sum(log(tmp.e2$values))
obj2 <- sum(log(tmp.e3$values)) + p
obj3 <- sum(diag(t3))
obj4 <- sum(diag(t5)) - 2 * sum(diag(t6))
objfun <- obj1 + obj2 + obj3 + obj4
loglik <- - (n/2) * objfun
if (asy == T) {
covMatrix <- tmp.xenv$covMatrix
asySE <- tmp.xenv$asySE
ratio <- tmp.xenv$ratio
}
}
else if (q == p) {
tmp.env <- env(X, Y, u, asy = TRUE, init = Pinit)
Gammahat <- tmp.env$Gamma
Gamma0hat <- tmp.env$Gamma0
Phihat <- diag(p)
Phi0hat <- NULL
etahat <- invsigX %*% crossprod(sigYX, Gammahat)
Omegahat <- tmp.env$Omega
Omega0hat <- tmp.env$Omega0
Deltahat <- sigX
Delta0hat <- NULL
betahat <- t(tmp.env$beta)
muhat <- colMeans(Y) - crossprod(betahat, colMeans(X))
SigmaYcXhat <- tmp.env$Sigma
SigmaXhat <- sigX
tmp.e1 <- eigen(Omegahat)
tmp.e2 <- eigen(Omega0hat)
tmp.e3 <- eigen(Deltahat)
invome <- chol2inv(chol(Omegahat))
invome0 <- chol2inv(chol(Omega0hat))
invdel <- chol2inv(chol(Deltahat))
t3 <- crossprod(Gammahat, sigY) %*% Gammahat %*% invome
t5 <- crossprod(etahat, sigX) %*% etahat %*% invome
t6 <- t(etahat) %*% crossprod(sigYX, Gammahat) %*% invome
obj1 <- sum(log(tmp.e1$values)) + sum(log(tmp.e2$values))
obj2 <- sum(log(tmp.e3$values)) + p
obj3 <- sum(diag(t3)) + (r - u)
obj4 <- sum(diag(t5)) - 2 * sum(diag(t6))
objfun <- obj1 + obj2 + obj3 + obj4
loglik <- - (n/2) * objfun
if (asy == T) {
M <- sigY - tcrossprod(betaOLS, sigYX)
covMatrix <- kronecker(M, invsigX)
asyFm <- matrix(sqrt(diag(covMatrix)), nrow = p)
m1 <- Gammahat %*% tcrossprod(Omegahat, Gammahat)
temp1 <- kronecker(m1, invsigX)
m2 <- crossprod(etahat, sigX) %*% etahat
m3 <- kronecker(invome0, m2) + kronecker(invome0, Omegahat)
m4 <- kronecker(Omega0hat, invome) - 2 * kronecker(diag(r - u), diag(u))
m5 <- m3 + m4
invm5 <- chol2inv(chol(m5))
m6 <- kronecker(Gamma0hat, etahat)
temp2 <- m6 %*% tcrossprod(invm5, m6)
covMatrix <- temp1 + temp2
asySE <- matrix(sqrt(diag(covMatrix)), nrow = p)
ratio <- asyFm/asySE
}
}
else{
tmp.stenv <- stenvMU(X, Y, q, u, initial1 = Pinit, initial2 = Ginit)
Gammahat <- tmp.stenv$Gammahat
Gamma0hat <- tmp.stenv$Gamma0hat
Phihat <- tmp.stenv$Phihat
Phi0hat <- tmp.stenv$Phi0hat
Deltahat <- crossprod(Phihat, sigX) %*% Phihat
Delta0hat <- crossprod(Phi0hat, sigX) %*% Phi0hat
invdel <- chol2inv(chol(Deltahat))
et1 <- t(Phihat) %*% crossprod(sigYX, Gammahat)
etahat <- invdel %*% et1
om1 <- sigYX %*% Phihat
om2 <- om1 %*% tcrossprod(invdel, om1)
om3 <- sigY - om2
Omegahat <- crossprod(Gammahat, om3) %*% Gammahat
Omega0hat <- crossprod(Gamma0hat, sigY) %*% Gamma0hat
betahat <- Phihat %*% tcrossprod(etahat, Gammahat)
muhat <- colMeans(Y) - crossprod(betahat, colMeans(X))
s1 <- Gammahat %*% tcrossprod(Omegahat, Gammahat)
SigmaYcXhat <- s1 + Gamma0hat %*% tcrossprod(Omega0hat, Gamma0hat)
s2 <- Phihat %*% tcrossprod(Deltahat, Phihat)
SigmaXhat <- s2 + Phi0hat %*% tcrossprod(Delta0hat, Phi0hat)
tmp.e1 <- eigen(Omegahat)
tmp.e2 <- eigen(Omega0hat)
tmp.e3 <- eigen(Deltahat)
tmp.e4 <- eigen(Delta0hat)
invome <- chol2inv(chol(Omegahat))
invome0 <- chol2inv(chol(Omega0hat))
invdel0 <- chol2inv(chol(Delta0hat))
t3 <- crossprod(Gammahat, sigY) %*% Gammahat %*% invome
peta <- Phihat %*% etahat
t5 <- crossprod(peta, sigX) %*% peta %*% invome
t6 <- t(peta) %*% crossprod(sigYX, Gammahat) %*% invome
obj1 <- sum(log(tmp.e1$values)) + sum(log(tmp.e2$values))
obj2 <- sum(log(tmp.e3$values)) + sum(log(tmp.e4$values)) + p
obj3 <- sum(diag(t3)) + (r - u)
obj4 <- sum(diag(t5)) - 2 * sum(diag(t6))
objfun <- obj1 + obj2 + obj3 + obj4
loglik <- - (n/2) * objfun
if (asy == T) {
M <- sigY - tcrossprod(betaOLS, sigYX)
covMatrix <- kronecker(M, invsigX)
asyFm <- matrix(sqrt(diag(covMatrix)), nrow = p)
m1 <- Gammahat %*% tcrossprod(Omegahat, Gammahat)
m2 <- Phihat %*% tcrossprod(invdel, Phihat)
temp1 <- kronecker(m1, m2)
d1 <- etahat %*% tcrossprod(invome, etahat)
m3 <- kronecker(d1, Delta0hat)
m4 <- kronecker(Deltahat, invdel0) + kronecker(invdel,
Delta0hat) - 2 * kronecker(diag(q), diag(p - q))
m5 <- m3 + m4
invm5 <- chol2inv(chol(m5))
m6 <- kronecker(tcrossprod(Gammahat, etahat), Phi0hat)
temp2 <- m6 %*% tcrossprod(invm5, m6)
d2 <- crossprod(etahat, Deltahat) %*% etahat
m8 <- kronecker(invome0, d2) + kronecker(invome0, Omegahat)
m9 <- kronecker(Omega0hat, invome) - 2 * kronecker(diag(r - u), diag(u))
m10 <- m8 + m9
invm10 <- chol2inv(chol(m10))
m11 <- kronecker(Gamma0hat, Phihat %*% etahat)
temp3 <- m11 %*% tcrossprod(invm10, m11)
covMatrix <- temp1 + temp2 + temp3
asySE <- matrix(sqrt(diag(covMatrix)), nrow = p)
ratio <- asyFm/asySE
}
}
return(list(beta = betahat, SigmaYcX = SigmaYcXhat, SigmaX = SigmaXhat,
Gamma = Gammahat, Gamma0 = Gamma0hat,
eta = etahat, Omega = Omegahat, Omega0 = Omega0hat,
mu = muhat, Phi = Phihat, Phi0 = Phi0hat,
Delta = Deltahat, Delta0 = Delta0hat,
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.