# 15 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 = length(lambdas), 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)
}