R/algo_leverage.R

In DavidJamesKent/stsci6520_hw2: STSCI6520 HW2

#' Leverage- or Uniform-subsample estimates for simple linear regression
#'
#' @param x Independent variable.
#' @param y Dependent variable.
#' @param r Size of subsample.
#' @param method Method, one of "uniform" or "leverage".
#' @return List containing estimates of intercept and slope, along with \code{samp}, the indices of the rows used in the subsample.
#' @export
#'
algo_leverage <- function(x,y,r,method="leverage") {

  n <- length(x)

  if(method=="leverage") {
    h <- (1/(sum(x^2))*(x^2))
    probs <- h/sum(h)
    samp <- sample(1:n,r,prob=h,replace=T)
    X <- cbind(rep(1,r),x[samp])
    Y <- y[samp]
    Phi_inv <- diag(1/probs[samp])
    B <- solve( t(X)%*% Phi_inv %*% X, t(X) %*% Phi_inv %*% Y )
    return(list(int=B[1],
                slope=B[2],
                samp=samp))
  } else if(method=="uniform") {
    probs <- rep(1/n,n)
    samp <- sample(1:n,r,prob=probs,replace=T)
    X <- cbind(rep(1,r),
               x[samp])
    Y <- y[samp]
    B <- solve(t(X)%*%X,t(X)%*%Y)
    return(list(int=B[1],
                slope=B[2],
                samp=samp))
  } else
    cat(paste("Method",method,"not recognized.\n"))
}