R/ls_solver.R

In andreabecsek/portfolio:

#' Least square solver
#'
#' @param X design matrix
#' @param y output vector
#'
#' @return estimated parameters
#'

rmse <- function(a,b) {sqrt(mean((a-b)^2))}

ridge <- function(X,y,lambda){
  n <- dim(X)[1]
  d <- dim(X)[2]
  XtX <- crossprod(X)
  w <- solve(XtX+lambda*diag(d),t(X)%*%y)
  pred <- X%*%w
  rmse <- rmse(pred,y)
  resid <- y - pred
  return(list(coef=w,pred=pred,resid=resid,rmse=rmse))
}


ls <- function(X,y){
  return(ridge(X,y,lambda=0))
}

# regression <- list(
#   ls = function(X,y) ls(X,y),
#   ridge = function(X,y,lambda) ridge(X,y,lambda))


# example
#generate data (multivariate)
set.seed(5)
n <- 100
x_1 <- runif(n, -10, 10)
x_2 <- runif(n, -10, 10)
X <- cbind(rep(1,n),x_1,x_2)
y <- 5*x_1-3*x_2 + rnorm(n)*4
## fit least squares model
fit <- ls(X,y)
plot(y,fit$resid)
plot(y,fit$pred)