R/algo_leverage.R

In sophiaycl/6520HW2: STSCI6520HW

algo_leverage = function(X, y, sampSize = floor(samplePct * length(y)),sampleType = 'Lev', samplePct = 0.1  ){
  n = length(y)

  LevPi = diag(tcrossprod(X, X))
  Phi = diag(rep(1, sampSize))

  UnifPi = rep(1, n)

  if (sampleType == 'Unif') {
    Sampled = sample(n, size = sampSize, prob = UnifPi)
  } else {
    Sampled = sample(n, size = sampSize, prob = LevPi)
    Phi = diag(1/LevPi[Sampled])
  }

  sampledX = X[Sampled,]
  sampledY = y[Sampled]

  betaHAT = crossprod(solve(crossprod(sampledX, Phi %*% sampledX)), crossprod(sampledX, Phi %*% as.matrix(sampledY)))

  return(betaHAT)
}

#' @param X The covariates matrix
#' @param y The response vector
#' @param sampSize  Number of subsamples
#' @param sampleType The type of subsampling rows
#' @param samplePct The percentage of subsamples
#'
#' @return The coefficient vector
#' @export
#'
#' @examples algo_leverage(matrix(rnorm(3000), nrow = 600), rnorm(600), sampleType = 'Unif')