40 lines
1.1 KiB
R
40 lines
1.1 KiB
R
#loocv
|
|
|
|
ridge.svd <- function(data, degre) {
|
|
X <- data$X
|
|
n <- nrow(X)
|
|
A <- scale(outer(c(X), 1:degre, "^"))
|
|
Ym <- mean(data$Y)
|
|
Y <- data$Y - Ym
|
|
As <- svd(A)
|
|
d <- As$d
|
|
function(lambda) {
|
|
coef <- c(As$v %*% ((d / (d^2 + lambda)) * (t(As$u) %*% Y)))
|
|
coef <- coef / attr(A,"scaled:scale")
|
|
inter <- Ym - coef %*% attr(A,"scaled:center")
|
|
coef <- c(inter, coef)
|
|
trace.H <- sum(d^2 / (d^2 + lambda))
|
|
pred <- polyeval(coef, X)
|
|
gcv <- sum( ((Y - pred) / (1 - (trace.H / n))) ^ 2 ) / n
|
|
list(coef = coef, gcv = gcv)
|
|
}
|
|
}
|
|
|
|
ridge.loocv <- function(data, deg, lambdas) {
|
|
errs <- double(length(lambdas))
|
|
coefs <- matrix(data = NA, nrow = nrow(data$X), ncol = 1+deg)
|
|
ridge <- ridge.svd(data, deg)
|
|
idx <- 1
|
|
for(lambda in lambdas) {
|
|
res <- ridge(lambda)
|
|
coefs[idx,] <- res$coef
|
|
errs[idx] <- res$gcv
|
|
idx <- idx + 1
|
|
}
|
|
err.min <- min(errs)
|
|
lambda.best <- lambdas[which(errs == err.min)]
|
|
coef.best <- coefs[which(errs == err.min),]
|
|
pred <- polyeval(coef.best, data$X)
|
|
mae <- mean(abs(pred - data$Y))
|
|
list(coef = coef.best, lambda = lambda.best, mae = mae)
|
|
} |