Nothing
lars.ss <-
function(y, X, ss, a.lst, S.v, C.v, col.plot, verbose=TRUE, plot=TRUE) {
check.ss.LARS(y, X, ss, a.lst, S.v, C.v, col.plot)
p <- dim(X)[2]
i <- 0
ind.v <- NULL
col.ind <- rep(1:p)
aic.v <- NULL
names.X <- colnames(X)
flag.stop <- FALSE
stack.ss <- NULL
A <- do.call(adiag, a.lst)
X <- as.matrix(X)
X.out <- as.matrix(X)
X.in <- matrix(0, nrow(X.out),1)
X.cand <- X%*%A
# Step 1: Initialize all regression coefficient estimates equal to 0, and let r=y-ybar.
if (verbose) {
cat("Step 1: Set r=y and initialize all betas to 0.", "\n")
}
r <- y # Note: ybar=0 since all variables were standardized to have mean=0 and SD=1.
beta.old <- array(0, dim=c(1,p), dimnames=list(NULL, names.X))
row.names(beta.old) <- i
# Step 2: Find the predictor among the p possible predictors that has the greatest absolute correlation with the residuals.
i <- i + 1
c.hat.v <- t(X.cand) %*%r
C.hat <- max(abs(c.hat.v))
index.maxabscor <- which.max(abs(c.hat.v))
ind.v <- col.ind[index.maxabscor]
if (verbose) {
cat("Step 2: Find the predictor that has the greatest absolute correlation with r.", "\n")
cat("\t", "Note that LARS picks the maximal absolute current correlation.", "\n", "\n")
}
a.lst.ind <- NULL
for (j in 1:length(a.lst)) {
if (length(which(unlist(dimnames(a.lst[[j]])[2])==colnames(X.out)[index.maxabscor]))==0) {
a.lst.ind[[j]] <- -1
} else {
a.lst.ind[[j]] <- which(unlist(dimnames(a.lst[[j]])[2])==colnames(X.out)[index.maxabscor])
}
}
index.lst <- which(a.lst.ind != -1)
index.col <- unlist(a.lst.ind[a.lst.ind != -1])
col.ind.v <- NULL
sub <- NULL
for (j in 1:S.v[index.lst]) {
col.ind.v[j] <- which(colnames(X)==dimnames(a.lst[[index.lst]])[2][[1]][j])
sub[j] <- which(col.ind==col.ind.v[j])
}
col.ind <- col.ind[-sub]
C.v[[index.lst]] <- index.col
a.lst[[index.lst]] <- a.lst[[index.lst]][,index.col, drop=FALSE]
if (S.v[index.lst]==1) {
stack.ss <- 1
} else {
stack.ss <- which(pmatch(ss, dimnames(a.lst[[index.lst]])[2])!="NA")
}
A <- do.call(adiag, a.lst)
X.cand <- X%*%A
if (verbose) {
cat("LARS Iteration ", i, sep="", "\n")
cat("\t", "Variable ", tail(ind.v, n=1), " added", sep="", "\n")
}
# Step 3: Compute c.hat.v and C.hat.
if (verbose) {
cat("\t", " Step 3: Compute c.hat.v and C.hat.", "\n")
}
c.hat.v <- t(X.cand) %*%r
print(c.hat.v)
C.hat <- max(abs(c.hat.v))
# Step 4: Set s.j=sign{c.hat.v}, and calculate X.in, A.in, u.in, and a.
if (verbose) {
cat("\t", " Set s.j=sign{c.hat.v}, and calculate X.in, A.in, u.in, and a.", "\n")
}
X.in <- X[,ind.v, drop=FALSE]
gamma.ind.v <- NULL
for (i in 1:dim(X.in)[2]) {
gamma.ind.v[i] <- which(colnames(X.cand)==colnames(X.in)[i])
}
X.in <- t(t(X.in)*sign(c.hat.v[gamma.ind.v]))
X.out <- X[,col.ind, drop=FALSE]
g.in <- t(X.in) %*%X.in
inv.g.in <- solve(t(X.in) %*%X.in)
A.in <- as.vector((t(rep(1,dim(X.in)[2])) %*% inv.g.in %*% rep(1,dim(X.in)[2]))^(-1/2))
w.in <- A.in * inv.g.in %*% rep(1,dim(X.in)[2])
u.in <- X.in %*% w.in
a <- t(X.cand) %*% u.in
# Step 5: Calculate gamma.hat and d.v.
if (verbose) {
cat("\t", " Step 5: Calculate gamma.hat and d.v.", "\n")
}
gamma.hat.v <- rep(0,p)
for (j in col.ind) {
comp1 <- (C.hat-(A%*%c.hat.v)[j])/(A.in-(A%*%a)[j])
comp2 <- (C.hat+(A%*%c.hat.v)[j])/(A.in+(A%*%a)[j])
comp <- c(comp1, comp2)
gamma.hat.v[j] <- min(comp[comp>0])
}
gamma.hat <- min(gamma.hat.v[gamma.hat.v>0])
gamma.hat.index <- which.min(gamma.hat.v[gamma.hat.v>0])
if (verbose) {
cat("\t", "\t", " Gamma values:", gamma.hat.v, "\n")
cat("\t", "\t", " Gamma hat:", gamma.hat, "\n")
}
sign.v <- sign(c.hat.v)
d.v <- rep(0,p)
d.v[ind.v] <- sign.v[gamma.ind.v] * w.in
# Step 6: Update beta.i <- beta.i-1 + gamma.hat*d.v, and update r.i <- y-X*beta.i.
beta.new <- beta.old[i,] + gamma.hat*d.v
if (verbose) {
cat("\t", " Step 6: Update betas and r.", "\n")
cat("\t", "\t", " Beta.new:", beta.new, "\n", "\n")
}
r <- y - X.cand%*%t(beta.new%*%A)
beta.old <- rbind(beta.old, beta.new)
beta.old
row.names(beta.old)[i+1] <- i
X.select <- X[,sort(ind.v)]
OLS <- lm(y~as.matrix(X.select)-1)
print(ind.v)
print(AIC(OLS))
aic.v[[i]] <- AIC(OLS)
while(flag.stop!=TRUE) {
i <- i + 1
# Repeat steps 3-6 until flag.stop=TRUE.
v1 <- pickvar.lars_v2.ss(r, y, X, X.cand, X.in, X.out, ss, a.lst, A, S.v, C.v, beta.old, ind.v, col.ind, gamma.hat.index, i, p, stack.ss, verbose)
r <- v1$r
X.in <- v1$X.in
X.out <- v1$X.out
X.cand <- v1$X.cand
a.lst <- v1$a.lst
A <- v1$A
C.v <- v1$C.v
beta.new <- v1$beta.new
beta.old <- rbind(beta.old, beta.new)
row.names(beta.old)[i+1] <- i
ind.v <- v1$ind.v
col.ind <- v1$col.ind
X.select <- X[,sort(ind.v)]
OLS <- lm(y~as.matrix(X.select)-1)
print(ind.v)
print(AIC(OLS))
aic.v[[i]] <- AIC(OLS)
gamma.hat.index <- v1$gamma.hat.index
stack.ss <- v1$stack.ss
flag.stop <- v1$flag.stop
if (flag.stop==TRUE) {
beta.nonzero <- beta.old%*%A
aic.sol <- which.min(aic.v)
beta.nonzero.aic <- beta.nonzero[(aic.sol+1),which(!beta.nonzero[(aic.sol+1),]==0)]
print(beta.nonzero[(i+1), ])
print(paste("No. of vars in model = ", dim(beta.nonzero)[2], sep=""))
print(paste("Total no. of vars considered = ", p, sep=""))
if (plot==TRUE) {
ind.stack <- cbind(ind.v, stack.ss)
stack.ss2 <- sortrows(ind.stack)[,2]
path.plot.lars_v2.ss(beta.nonzero, ind.v, i, stack.ss2, aic.sol, col.plot)
}
break
}
}
return(list(beta=beta.nonzero, beta.aic=beta.nonzero.aic, ind.v=ind.v, aic.v=aic.v, stack.ss=stack.ss))
}
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.