From 3dafda091b49fade2c5a40a6db12485243337ff1 Mon Sep 17 00:00:00 2001
From: Pierre-Edouard Portier <p6e7p7@freeshell.org>
Date: Thu, 2 Dec 2021 18:43:18 +0100
Subject: [PATCH] =?UTF-8?q?Impl=C3=A9mentation=20de=20la=20validation=20cr?=
 =?UTF-8?q?ois=C3=A9e=20un=20contre=20tous=20(loocv)=20pour=20un=20mod?=
 =?UTF-8?q?=C3=A8le=20de=20r=C3=A9gression=20ridge.?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 15_loocv.R   | 39 +++++++++++++++++++++++++++++++++++++++
 15_loocv.Rmd | 27 +++++++++++++++++++++++++++
 2 files changed, 66 insertions(+)

diff --git a/15_loocv.R b/15_loocv.R
index 798ce17..8e7fd3c 100644
--- a/15_loocv.R
+++ b/15_loocv.R
@@ -1 +1,40 @@
 #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)
+}
\ No newline at end of file
diff --git a/15_loocv.Rmd b/15_loocv.Rmd
index 0ab9b1c..3517c9a 100644
--- a/15_loocv.Rmd
+++ b/15_loocv.Rmd
@@ -134,3 +134,30 @@ 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
+deg <- 8
+data = gendat(n,0.2)
+splitres <- splitdata(data,0.8)
+entr <- splitres$entr
+test <- splitres$test
+lambdas <- c(1E-8, 1E-7, 1E-6, 1E-5, 1E-4, 1E-3, 1E-2, 1E-1, 1)
+ridge.res <- ridge.loocv(entr, deg, lambdas)
+plt(entr,f)
+pltpoly(ridge.res$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 contre tous un polynôme de degré au plus égal à `r deg` qui modélise au mieux ces données. La valeur de $\lambda$ retenue est : `r ridge.res$lambda` et l'erreur absolue moyenne sur le jeu d'entraînement est : `r ridge.res$mae`.
+
+```{r}
+testpred <- polyeval(ridge.res$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(ridge.res$coef)
+```