R/elnet_coord.R

In sophiaycl/6520HW2: STSCI6520HW

#' Elastic net regression
#' @examples elnet_coord(matrix(rnorm(400), nrow = 100), runif(100), beta0 = rep(1,4))
elnet_coord = function(X, Y,lambda = 1, a = 0.5,  iterMax = 100, tol = 10^(-8), GLMNET = FALSE, beta0 = rep(0, dim(X)[2])) {
  if (!GLMNET) {
    n = length(Y)
    for (i in 1:iterMax) {
      updateBeta = beta0
      for (j in 1:dim(X)[2]) {
        rho = 2 * crossprod(X[,j], (Y - crossprod(t(X[,-j]), updateBeta[-j])))/n
        z = 2 * crossprod(X[j,], X[j,])/n + 2 * lambda * (1 - a)
        if (rho < - a * lambda) {
          updateBeta[j] = (rho + a * lambda)/z
        } else {
          if (rho > a * lambda) {
            updateBeta[j] = (rho - a * lambda)/z
          } else {updateBeta[j] = 0}
        }
      }
      if (norm(as.matrix(updateBeta - beta0)) < tol) {
        break
      }
      beta0 = updateBeta
    }
    return(updateBeta)
  } else {
    return(glmnet::glmnet(X, Y, alpha = a, lambda = lambda, thresh = tol, maxit = iterMax))
  }
}
#' @param X Covariate matrix
#' @param Y Response vector
#' @param GLMNET An indicator whether

#' @return If \code{GLMNET = FALSE}
#' @export
##' @param a
#' @param lambda  parameter of the weight magnitude of penalty terms.
#' @param tol
#' @param iterMax Maximum iteration times
#' @param beta0 Initial guess for coefficients

#' @seealso \code{\link[glmnet]{glmnet}}