Nothing
lad <-
function(X, y, numdir=NULL, nslices=NULL, numdir.test=FALSE,...)
{
mf <- match.call();
op <- dim(X); p <- op[2]; n <- op[1]; yy <- y; yfactor <- FALSE
verbose <- mf$verbose; if (is.null(verbose)) verbose=FALSE;
if (is.factor(y)) {yy <- as.vector(as.integer(factor(y, levels=unique(y)))); yfactor <- TRUE}
if (n != length(yy)) stop("X and y do not have the same number of observations")
if (p > n) stop("This method requires n larger than p");
vnames <- dimnames(X)[[2]];
if (is.null(vnames)) vnames <- paste("X", 1:p, sep="");
X <- as.matrix(X);
Sigmatilde <- cov(X);
if (is.null(nslices)) if (verbose) message("The response is treated as categorical.")
if (is.null(nslices))
{
groups <- unique(sort(yy));
use.nslices <- length(groups);
if (identical(use.nslices, 1)) stop("The response has a single level.")
numdir <- min(p, use.nslices-1);
N_y <-vector(length=use.nslices);
Deltatilde_y = vector(length=use.nslices, "list");
for (i in 1:use.nslices)
{
N_y[i] <- length(yy[yy==groups[i]]);
X_y <- X[which(yy == groups[i]),];
Deltatilde_y[[i]] <- cov(X_y);
}
}
if (!is.null(nslices))
{
if (identical(nslices, 1)) stop("You need at least two slices.")
use.nslices=nslices;
if (is.null(numdir)) numdir <- min(p, nslices-1);
ysort <- ldr.slices(yy, nslices);
N_y <-vector(length=nslices);
Deltatilde_y <- vector(length=nslices, "list");
for (i in 1:nslices)
{
N_y[i] <- sum(ysort$slice.indicator==i);
X_y <- X[which(ysort$slice.indicator == i),];
Deltatilde_y[[i]] <- cov(X_y);
}
}
Sigmas <- list(Sigmatilde=Sigmatilde, Deltatilde_y=Deltatilde_y, N_y=N_y);
Deltatilde <- 0;
for (i in 1:use.nslices) Deltatilde <- Deltatilde + N_y[i] * Deltatilde_y[[i]];
InitialMatrix <- function(X, y, nslices)
{
op <- dim(X); n <- op[1]; p <- op[2];
if (is.null(nslices))
{
yunique <- sort(unique(y));
ygroups <- y;
nslices=length(unique(y));
}
else
{
Slicing <- ldr.slices(y, nslices=nslices);
yunique <- sort(unique(Slicing$slice.indicator));
ygroups <- Slicing$slice.indicator;
}
SigmaX <- cov(X);
X0 <- as.matrix(X - kronecker(cbind(rep(1,n)), rbind(apply(X,2,mean))));
svdS<-svd(SigmaX);
InvSqrtX <- svdS$u %*% diag(c(svdS$d)^(-0.5)) %*% t(svdS$v);
Z <- X0%*%InvSqrtX;
pv <- vector(length=nslices);
for (i in 1:nslices) pv[i] <- sum(ygroups==yunique[i])/length(y);
cez <- matrix(0,p,nslices);
cezz <- array(0,dim=c(p,p,nslices));
for (i in 1:nslices){
for (j in 1:n){
cez[,i] <- cez[,i] + Z[j,]*(ygroups[j]==i);
cezz[,,i] <- cezz[,,i] + Z[j,]%*%t(Z[j,])*(ygroups[j]==i);
}
cez[,i] <- cez[,i]/n/pv[i];
cezz[,,i] <- cezz[,,i]/n/pv[i];
}
mat1 <- mat2 <- mat4 <- mat5 <- matrix(0,p,p); mat3 <- 0;
for (i in 1:nslices){
mat1 <- mat1 + cezz[,,i]%*%cezz[,,i]*pv[i];
mat2 <- mat2 + cez[,i]%*%t(cez[,i])*pv[i];
mat3 <- mat3 + sum(cez[,i]^2)*pv[i];
mat4 <- mat4 + cezz[,,i]%*%cez[,i]%*%t(cez[,i])*pv[i];
mat5 <- mat5 + sum(cez[,i]^2)*cez[,i]%*%t(cez[,i])*pv[i];
}
dr <- 2*mat1 + 2*mat2%*%mat2 + 2*mat3*mat2 - 2*diag(p);
vector.dr <- eigen(dr)$vectors;
dr.direction <- InvSqrtX %*%vector.dr;
Wdr <- orthonorm(dr.direction);
# using SIR
sir <- mat2; vector.sir <- eigen(sir)$vectors;
sir.direction <- InvSqrtX %*%vector.sir;
Wsir <- orthonorm(sir.direction);
# Using SAVE
save <- mat1 - diag(p) - mat4 - t(mat4) + 2*sir + mat5;
vector.save <- eigen(save)$vectors;
save.direction <- InvSqrtX %*%vector.save;
Wsave <- orthonorm(save.direction);
return(list(Wsir=Wsir, Wsave=Wsave, Wdr=Wdr))
}
onelad <-function(d,...)
{
if (d==0)
{
Gammahat <- NULL;
Deltahat <- Deltatilde;
terme0 <- -(n*p)*(1+log(2*pi))/2;
terme1 <- -n*log( det(Sigmatilde) )/2;
loglik = terme0 + terme1;
Deltahat_y <- Deltatilde_y;
numpar <- p + p*(p+1)/2;
AIC <- -2*loglik + 2 * numpar;
BIC <- -2*loglik + log(n) * numpar;
}
else if (d==p)
{
Gammahat=diag(p);
Deltahat=Deltatilde;
terme0 <- -(n*p)*(1+log(2*pi))/2;
terme3<-0;
for (i in 1:length(N_y))
{
terme3<- terme3 - N_y[i] * log( det(Deltatilde_y[[i]]) )/2;
}
loglik = terme0 + terme3;
Deltahat_y <- Deltatilde_y;
numpar <- p+(use.nslices-1)*d + p*(p+1)/2 + (use.nslices-1)*d*(d+1)/2;
AIC <- -2*loglik + 2 * numpar;
BIC <- -2*loglik + log(n) * numpar;
}
else
{
objfun <- function(W)
{
Q <- W$Qt; d <- W$dim[1]; p <- W$dim[2];
Sigmas <- W$Sigmas;
n <- sum(Sigmas$N_y);
U <- matrix(Q[,1:d], ncol=d);
V <- matrix(Q[,(d+1):p], ncol=(p-d));
# Objective function
terme0 <- -(n*p)*(1+log(2*pi))/2;
terme1 <- -n*log(det(Sigmas$Sigmatilde))/2;
term <- det(t(U)%*%Sigmas$Sigmatilde%*%U);
if (is.null(term)|is.na(term)|(term <= 0)) return(NA) else terme2 <- n*log(term)/2 ;
terme3<-0;
for (i in 1:length(Sigmas$N_y))
{
term <- det(t(U)%*%Sigmas$Deltatilde_y[[i]]%*%U);
if (is.null(term)|is.na(term)|(term <=0)) {break; return(NA)};
terme3 <- terme3 - Sigmas$N_y[i] * log(term)/2;
}
value <- terme0 + terme1 + terme2 + terme3;
#Gradient
terme0 <- solve(t(U)%*%Sigmas$Sigmatilde%*%U);
terme1 <- sum(Sigmas$N_y) * t(V)%*%Sigmas$Sigmatilde%*%U%*%terme0;
terme2 <- 0;
for (i in 1:length(Sigmas$N_y))
{
temp0 <- solve(t(U)%*%Sigmas$Deltatilde_y[[i]]%*%U);
temp <-Sigmas$N_y[i]*t(V)%*%Sigmas$Deltatilde_y[[i]]%*%U%*%temp0;
terme2 <- terme2 - temp;
}
gradient <- t(terme1 + terme2);
return(list(value=value, gradient=gradient))
}
objfun <- assign("objfun", objfun, envir=.BaseNamespaceEnv);
if (is.null(mf$sim_anneal))
{
Wlist <- InitialMatrix(X=X, y=yy, use.nslices);
values <- c(objfun(W=list(Qt=Wlist$Wsir, dim=c(d,p), Sigmas=Sigmas))$value,
objfun(W=list(Qt=Wlist$Wsave, dim=c(d,p), Sigmas=Sigmas))$value,
objfun(W=list(Qt=Wlist$Wdr, dim=c(d,p), Sigmas=Sigmas))$value);
W <- list(Qt=Wlist[[which(values == max(values))[1]]], dim=c(d,p), Sigmas=Sigmas)
grassmann <- GrassmannOptim(objfun, W, ...);
}
if (!is.null(mf$sim_anneal))
{
W <- list(dim=c(d,p), Sigmas=Sigmas);
grassmann <- GrassmannOptim(objfun, W,...);
}
Gammahat = matrix(grassmann$Qt[,1:d], ncol=d);
dimnames(Gammahat) <- list(vnames, paste("Dir", 1:d, sep=""))
invDeltahat <- solve(Sigmatilde) +
Gammahat%*%solve(t(Gammahat)%*%Deltatilde%*%Gammahat)%*%t(Gammahat) -
Gammahat%*%solve(t(Gammahat)%*%Sigmatilde%*%Gammahat)%*%t(Gammahat);
Deltahat <- solve(invDeltahat); Pj <- projection(Gammahat, Deltahat);
Deltahat_y <- vector("list", use.nslices);
for (i in 1:use.nslices)
{
Deltahat_y[[i]] <- Deltahat + t(Pj) %*% (Deltatilde_y[[i]] - Deltahat)%*%Pj;
}
terme0 <- -(n*p)*(1+log(2*pi))/2;
terme1 <- -n*log( det(Sigmatilde) )/2;
terme2 <- n*log(det(t(Gammahat)%*%Sigmatilde%*%Gammahat))/2;
terme3<-0;
for (i in 1:length(N_y))
{
terme3 <- terme3 - N_y[i] * log( det(t(Gammahat)%*%Deltatilde_y[[i]]%*%Gammahat) )/2;
}
loglik = terme0 + terme1 + terme2 + terme3;
numpar <- p+(use.nslices-1)*d+p*(p+1)/2+d*(p-d)+(use.nslices-1)*d*(d+1)/2;
AIC <- -2*loglik + 2 * numpar;
BIC <- -2*loglik + log(n) * numpar;
}
return(list(Gammahat=Gammahat, Deltahat=Deltahat, Deltahat_y=Deltahat_y, loglik=loglik,
numpar=numpar, aic=AIC, bic=BIC))
}
if (identical(numdir.test, FALSE))
{
fit <- onelad(numdir, ...);
R <- X%*%fit$Gammahat;
ans <- list(R=R, Gammahat=fit$Gammahat, Deltahat=fit$Deltahat, Deltahat_y=fit$Deltahat_y,
loglik=fit$loglik, aic=fit$aic, bic=fit$bic, numpar=fit$numpar, numdir=numdir,
yfactor=yfactor, y=y, model="lad", call=match.call(expand.dots=TRUE),numdir.test=numdir.test);
class(ans) <- "lad";
return(invisible(ans));
}
aic <- bic <- numpar <- loglik <- vector(length=numdir+1);
Deltahat <- Deltahat_y <- Gammahat <-vector("list");
for (i in 0:numdir)
{
fit <- onelad(i, ...);
Gammahat[[i+1]] <-fit$Gammahat;
Deltahat[[i+1]] <- fit$Deltahat;
Deltahat_y[[i+1]] <- fit$Deltahat_y;
loglik[i+1] <- fit$loglik;
numpar[i+1] <- fit$numpar;
aic[i+1] <- fit$aic;
bic[i+1] <- fit$bic;
}
R <- X%*%Gammahat[[numdir+1]];
ans <- list(R=R, Gammahat=Gammahat, Deltahat=Deltahat, Deltahat_y=Deltahat_y, loglik=loglik, aic=aic,
bic=bic, numpar=numpar, numdir=numdir, model="lad", call=match.call(expand.dots=TRUE),
yfactor=yfactor, y=y, numdir.test=numdir.test);
class(ans)<- "lad";
return(invisible(ans))
}
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.