Nothing
senvMU <- function(X, Y, u, initial = 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.")
objfun <- function(d, Gamma, X, Y){
X <- as.matrix(X)
Y <- as.matrix(Y)
a <- dim(Y)
n <- a[1]
r <- a[2]
p <- ncol(X)
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))
invsigY <- chol2inv(chol(sigY))
betaOLS <- sigYX %*% invsigX
U <- tcrossprod(betaOLS, sigYX)
M <- sigY - U
d1 <- c(1, d)
Lambda <- diag(d1)
invLambda <- diag(1 / d1)
m1 <- crossprod(Gamma, Lambda)
eigtem1 <- eigen(m1 %*% tcrossprod(invsigY, m1))
m2 <- crossprod(Gamma, invLambda)
eigtem2 <- eigen(m2 %*% tcrossprod(M, m2))
temp1 <- sum(log(eigtem1$values))
temp2 <- sum(log(eigtem2$values))
objfun <- temp1 + temp2
return(objfun)
}
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))
invsigY <- chol2inv(chol(sigY))
betaOLS <- sigYX %*% invsigX
U <- tcrossprod(betaOLS, sigYX)
M <- sigY - U
tmp <- envMU(M, U, u, initial = initial)
Gamma <- tmp$Gammahat
d.init <- rep(1, (r - 1))
k2 <- rep(0, (r - 1))
tmp.init <- Rsolnp::solnp(pars = d.init, fun = objfun, LB = k2,
control = list(delta = 1e-8, tol = 1e-6, trace = 0),
Gamma = Gamma, X = X, Y = Y)
d <- tmp.init$pars
obj1 <- objfun(d, Gamma, X, Y)
maxiter = 20000
ftol = 1e-04
i <- 1
while (i < maxiter) {
d1 <- c(1, d)
L1 <- diag(d1)
invL1 <- diag(1 / d1)
M1 <- invL1 %*% M %*% invL1
U1 <- invL1 %*% U %*% invL1
tmp1 <- envMU(M1, U1, u)
Gamma <- tmp1$Gammahat
temp1 <- Rsolnp::solnp(pars = d, fun = objfun, LB = k2,
control = list(delta = 1e-8, tol = 1e-6, trace = 0),
Gamma = Gamma, X = X, Y = Y)
d <- temp1$pars
obj2 <- objfun(d, Gamma, X, Y)
if (abs(obj1 - obj2) < ftol) {
break
}
else {
obj1 <- obj2
i <- i + 1
}
}
Gamma0 <- tmp1$Gamma0hat
d1 <- c(1, d)
Lambda <- diag(d1)
return(list(Gammahat = Gamma, Gamma0hat = Gamma0, Lambdahat = Lambda,
objfun = obj2))
}
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.