Skip to content

Instantly share code, notes, and snippets.

@MGallow MGallow/linear.R
Created May 19, 2017

Embed
What would you like to do?
weighted ridge regression script
## weighted ridge regression
linearr = function(X, y, lam = 0, weights = NULL, intercept = TRUE,
kernel = FALSE) {
# checks
n = dim(X)[1]
p = dim(X)[2]
if (is.null(weights)) {
weights = rep(1, n)
}
if (length(weights) != n)
stop("weights must be length ", n)
if (length(lam) > 1)
stop("lam must be a scalar!")
if (lam < 0)
stop("lam must be nonnegative!")
if ((kernel == TRUE) & (lam == 0))
stop("must specify lam to use kernel!")
# initialization
X. = as.matrix(X)
y. = as.matrix(y)
if (intercept == TRUE) {
# if first column of ones, then remove it
if (all(X.[, 1] == rep(1, n))) {
X. = X.[, -1]
}
# center the data
X_bar = (as.matrix(t(weights)) %*% X.)/sum(weights)
y_bar = sum(weights * y.)/sum(weights)
X. = X. - rep(1, n) %*% X_bar
y. = y. - y_bar
}
root = sqrt(weights)
X. = root * X.
y. = root * y.
# SVD
svd = svd(X.)
# adjust d vector for regularization and diagonalize
d_adj = diag(svd$d/(svd$d^2 + lam))
# calculate beta estimates
betas = svd$v %*% d_adj %*% t(svd$u) %*% y.
rownames(betas) = NULL
# add intercept, if needed
if (intercept == TRUE) {
# calculate intercept
b1 = y_bar - X_bar %*% betas
rownames(b1) = c("intercept")
betas = rbind(b1, betas)
}
returns = list(coefficients = betas)
return(returns)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.