hwk 6 no. 4

gistfile1.txt

data(prostate)
prostate <- prostate[,1:9]

#scale all the predictors
data.st <- scale(prostate[,-9])
data.st <- sapply(1:8, function(x) data.st[,x]/fxn.4$norm(data.st[,x]))
data.st <- cbind(data.st, prostate[,9])

#roughly unit vectors
sapply(1:9, function(x) fxn.4$norm(data.st[,x]))

#roughly mean of 0
colMeans(data.st)

#so, per the book, this means alpha.0 = mean(y)
alphas <- rep(0, 9)
alphas[1] <- mean(data.st[,9])

fxn.4 <- list(
  norm = function(x) return(sqrt(sum(x)^2)),
  thresh =function(z, lambda) {
  sign <- ifelse(z <= 0, -1, 1)
  t <- abs(z) - lambda
  step <- ifelse( t > 0, t, 0)
  return(sign*step)
  }
  )


# coordinate descent algorithm
c.descent <- function(start = rep(0,8), fxn = fxn.4, y.i = data.st[,9], x.i = data.st[,1:8], lambda = .1, eps = .01) {

# start out with the initial parameter values
point <- start
alphas <- rep(NULL, length(point))

while(fxn$norm(alphas - point) > eps) {
point <- alphas

# iterate through the predictors
for(i in 1:dim(x.i)[2]) {
# calculate the residuals based on all other predictors 
r <- y.i - x.i[,-i] %*% point[-i]

# calculate the least squares coefficient of the residuals regressed on target predictor
alpha.t <- coefficients(lm(r ~ x.i[,i] - 1))[1]

# soft thresholding step
# store the values in the working point and the next point
alphas[i] <- fxn$thresh(alpha.t, lambda)
}
}
return(alphas) 
}