K-fold-cv pour la régression ridge basée sur le svd de la matrice des données.

14_geometrie_ridge_svd.R

@@ -25,9 +25,15 @@ ridge.gram <- function(alpha, data, degre) {
   coef <- c(inter, coef)
 }
 
-ridge.svd <- function(data, degre) {
-  A <- scale(outer(c(data$X), 1:degre, "^"))
+ridge.svd <- function(data, degre, fold = FALSE) {
+  if (length(fold) == 1 && fold == FALSE) {
+    X <- data$X
     Y <- data$Y
+  } else {
+    X <- data$X[!fold,]
+    Y <- data$Y[!fold]
+  }
+  A <- scale(outer(c(X), 1:degre, "^"))
   Ym <- mean(Y)
   Y <- Y - Ym
   As <- svd(A)
@@ -58,20 +64,23 @@ kfoldridge <- function(K, alphas, data, degre) {
   folds <- sample(folds, N)
   maes <- matrix(data = NA, nrow = K, ncol = length(alphas))
   colnames(maes) <- alphas
-  alpha_idx <- 1
-  for(alpha in alphas) {
   for(k in 1:K) {
     fold <- folds == k
-      coef <- ridge(alpha, data, degre, fold)
-      pred <- polyeval(coef, data$X[fold,])
-      maes[k,alpha_idx] <- mean(abs(pred - data$Y[fold]))
-    }
+    ridge <- ridge.svd(data, degre, fold)
+    alpha_idx <- 1
+    X <- data$X[fold,]
+    Y <- data$Y[fold]
+    for(alpha in alphas) {
+      coef <- ridge(alpha)
+      pred <- polyeval(coef, X)
+      maes[k,alpha_idx] <- mean(abs(pred - Y))
       alpha_idx <- alpha_idx + 1
     }
+  }
   mmaes <- colMeans(maes)
   minmmaes <- min(mmaes)
   bestalpha <- alphas[which(mmaes == minmmaes)]
-  fold <- folds == K+1 # vector of FALSE
-  coef <- ridge(bestalpha, data, degre, fold)
+  ridge <- ridge.svd(data, degre)
+  coef <- ridge(bestalpha)
   list(coef = coef, maes = maes, alpha = bestalpha)
 }

14_geometrie_ridge_svd.Rmd

@@ -73,7 +73,34 @@ Nous vérifions que les coefficients trouvés par les deux approches sont identi
 all.equal(coef.gram,coef.svd)
 ```
 
-La décomposition en valeurs singulières de $\mathbf{X}$ donne les coefficients de la régression Ridge pour toutes les valeurs possibles du coefficient de régularisation $\lambda$.
+La décomposition en valeurs singulières de $\mathbf{X}$ donne les coefficients de la régression Ridge pour toutes les valeurs possibles du coefficient de régularisation $\lambda$. Nous pouvons donc implémenter de façon efficace une validation croisée à k plis.
+
+```{r}
+alphas <- c(1E-8, 1E-7, 1E-6, 1E-5, 1E-4, 1E-3, 1E-2, 1E-1, 1)
+reskfold <- kfoldridge(K = 10, alphas = alphas, data = entr, degre = deg1)
+plt(entr,f)
+pltpoly(reskfold$coef)
+```
+
+Ci-dessus, nous avons généré un jeu de données composé de `r n` observations et nous avons calculé par validation croisée un polynôme de degré au plus égal à `r deg1` qui modélise au mieux ces données. La valeur de $\alpha$ retenue est : `r reskfold$alpha`.
+
+Traçons un boxplot des erreurs commises sur les plis de validation pour chaque valeur de $\alpha$.
+
+```{r}
+boxplot(reskfold$maes)
+```
+
+```{r}
+testpred <- polyeval(reskfold$coef, test$X)
+testmae <- mean(abs(testpred - test$Y))
+```
+
+Ce meilleur modèle atteint une erreur absolue moyenne de `r testmae` sur le jeu de test.
+
+```{r}
+plt(test,f)
+pltpoly(reskfold$coef)
+```
 
 # Régression ridge et géométrie
 

15_loocv.R

@@ -1,30 +1 @@
 #loocv
-
-# Résolution d'un système linéaire correspondant à la matrice de Gram pour
-# un polynôme de degré fixé et avec l'ajout d'un facteur de régularisation en
-# norme L2 dont l'importance est contrôlée par l'hyperparamètre alpha.
-ridge.gram <- function(alpha, data, degre) {
-  A <- scale(outer(c(data$X), 1:degre, "^"))
-  Y <- data$Y
-  Ym <- mean(Y)
-  Y <- Y - Ym
-  gram <- t(A) %*% A
-  diag(gram) <- diag(gram) + alpha
-  coef <- solve(gram, as.vector(t(A) %*% Y))
-  coef <- coef / attr(A,"scaled:scale")
-  inter <- Ym - coef %*% attr(A,"scaled:center")
-  coef <- c(inter, coef)
-}
-
-ridge.svd <- function(alpha, data, degre) {
-  A <- scale(outer(c(data$X), 1:degre, "^"))
-  Y <- data$Y
-  Ym <- mean(Y)
-  Y <- Y - Ym
-  As <- svd(A)
-  d <- As$d
-  coef <- c(As$v %*% ((d / (d^2 + alpha)) * (t(As$u) %*% Y)))
-  coef <- coef / attr(A,"scaled:scale")
-  inter <- Ym - coef %*% attr(A,"scaled:center")
-  coef <- c(inter, coef)
-}

15_loocv.Rmd

@@ -134,23 +134,3 @@ Nous remarquons que cette mesure de l'erreur peut être instable quand au moins
  &\text{avec } tr(\mathbf{H}(\lambda)) = \sum_{d_j>0} \frac{d_j^2}{d_j^2 + \lambda}
 \end{align*}
 
-```{r}
-set.seed(1123)
-n <- 100
-deg1 <- 8
-data = gendat(n,0.2)
-splitres <- splitdata(data,0.8)
-entr <- splitres$entr
-test <- splitres$test
-alpha <- 1E-5
-coef.gram <- ridge.gram(alpha, entr, deg1)
-plt(entr,f)
-pltpoly(coef.gram)
-```
-
-```{r}
-coef.svd <- ridge.svd(alpha, entr, deg1)
-plt(entr,f)
-pltpoly(coef.svd)
-```
-