R/SLOB.R

In StatsIMUWr/slobeC: What the package does (short line)

#' Find SLOB estimator
#'
#' @param start initial beta vector
#' @param Y response variable
#' @param X data matrix
#' @param lambda0 initial lambda vector
#' @param a ?
#' @param b ?
#' @param max_iter Maximum number of iterations
#' @param tolerance Tolerance
#' 
#' @author Malgorzata Bogdan
#'
#' @return list that consists of...
#'
#' @export
#'

SLOB <- function(start, Y, X, lambda0, a, b, max_iter = 100, tolerance = 1e-4) {
  eps<-1
  p <- ncol(X)
  n <- nrow(X)

  beta_k <- start
  beta_new <- beta_k

  o <- p - rank(abs(beta_k), ties.method = "max") + 1
  sigma <- sd(Y)
  lambdan <- lambda0[o] * sigma
  lambdan_sigma_inv <- lambda0[o] / sigma
  lambda_sigma <- lambda0 * sigma
  c_k <- min((sum(abs(start) > 0) + 1) / sum(abs(start[start != 0])) / lambda_sigma[p], 1)
  if(!is.finite(c_k)) {
    c_k <- 1
  }
  theta <- (sum(start != 0) + a) / (p + b + a)
  pstar_k <- 1 / (1 + (1 - theta) / theta /c_k * exp(-abs(beta_k)*lambdan_sigma_inv * (1 - c_k)))

  k <- 1
  while(eps > tolerance & k < max_iter) {
    wts <- pstar_k * c_k + (1 - pstar_k)
    revwts <- 1 / wts
    revwts <- as.vector(revwts)
    Xtemp <- sweep(X, 2, revwts, '*')
    utemp <- slope_admm(Xtemp, Y, rep(0, p), rep(0, p),
                        lambda_seq = lambda_sigma, 1)
    utemp1 <- utemp[["z"]]
    beta_k <- revwts * utemp1
    o1 <- p- rank(abs(utemp1), ties.method = "max") + 1
    RSS <- sum((Y - X %*% beta_k) ^ 2)
    sum_lamwbeta <- sum(lambda0[o1] * abs(utemp1))
    sigma <- (sqrt(sum_lamwbeta ^ 2 + 4 * n * RSS) + sum_lamwbeta) / (2 * n)
    lambda_sigma <- lambda0 * sigma
    lambda_sigma_inv <- lambda0 / sigma
    pstar_k <- 1 / (1 + (1 - theta) / theta / c_k * exp(-abs(beta_k)*lambda_sigma_inv[o1]*(1 - c_k)))
    rate_temp <- sum(abs(beta_k) * lambda_sigma_inv[o1] * pstar_k)
    c_k <- min((sum(pstar_k) + 1) / rate_temp, 1)
    if(!is.finite(c_k)) {
      c_k<-1
    }
    theta <- (sum(pstar_k) + a)/(p + b + a)
    k <- k + 1
    eps <- sum((beta_new - beta_k)^2)
    beta_new <- beta_k
  }

  list(beta = beta_k, pgamma = pstar_k, sigma = sigma, c = c_k, theta = theta)
}