Progrès...

ML1_regression.Rmd

@@ -1,6 +1,9 @@
 ---
 title: "Régression"
-output: html_notebook
+output:
+  html_notebook:
+    number_sections: yes
+  pdf_document: default
 ---
 
 \renewcommand{\vec}[1]{\mathbf{#1}}
@@ -79,33 +82,39 @@ gendat <- function(N, sd) {
   # sd: écart type d'un bruit gaussien de moyenne nulle
   X = runif(N)
   Y = f(X) + rnorm(N, mean=0, sd=sd)
+  dim(X) <- c(N,1) # en général chaque observation est décrite par plusieurs variables
+                   # et X est une matrice avec autant de lignes que d'observations
+                   # et autant de colonnes que de variables. Sur notre exemple,
+                   # chaque observation n'est décrite que par une seule variable.
   list(X = X, Y = Y)
 }
 ```
 On se donne une fonction pour afficher simultanément un jeu de données et la fonction utilisée pour le générer.
 ```{r}
-plt <- function(data, f) {
+plt <- function(data, f, ...) {
   xs = seq(0,1,length.out=100)
-  plot(xs, f(xs), type="l")
+  plot(xs, f(xs), type="l", ...)
   points(data$X, data$Y)
 }
 ```
 On affiche un exemple de jeu de données.
 ```{r}
 set.seed(1123)
-data = gendat(5,0.2)
+data = gendat(10,0.2)
 plt(data,f)
 ```
-La fonction ci-dessous résout le système linéaire correpsondant à la matrice de Vandermonde et retourne les coefficients du polynôme résultat.
+La fonction ci-dessous résout le système linéaire correspondant à la matrice de Vandermonde et retourne les coefficients du polynôme résultat.
 ```{r}
 polyreg1 <- function(data) {
-  vandermonde = outer(data$X, 0:(length(data$X)-1), "^")
+  xs <- c(data$X) # on transforme la matrice X de dimension Nx1 sur notre exemple,
+                  # en un vecteur
+  vandermonde = outer(xs, 0:(length(xs)-1), "^")
   solve(vandermonde, data$Y)
 }
 ```
 La fonction ci-dessous évalue un polynôme en un point `x` étant donné ses coefficients `coef`.
 ```{r}
-polyeval <- function(x,coef) {
+polyeval <- function(coef,x) {
   powers = 0:(length(coef)-1)
   f = function(y) { sum(coef * y^powers) }
   sapply(x,f)
@@ -115,7 +124,7 @@ La fonction ci-dessous ajoute au graphe courant le tracé en pointillés d'un po
 ```{r}
 pltpoly <- function(coef) {
   xs = seq(0,1,length.out=100)
-  lines(xs, polyeval(xs,coef), lty="dotted")
+  lines(xs, polyeval(coef,xs), lty="dotted")
 }
 ```
 On affiche la fonction génératrice, le jeu de donnée et le polynôme qui passe par chaque point du jeu de données.
@@ -169,6 +178,7 @@ $$
  & \vec{x}^T\vec{A}^T\vec{A}\vec{x} - 2\vec{x}^T\vec{A}^T\vec{b} + \vec{b}^T\vec{b}
 \end{align*}
 $$
+
 Cette dernière expression quadratique en $\vec{x}$ correspond à une surface convexe. Donc, le minimum de cette expression peut être calculé en annulant sa dérivée (penser à une courbe $y = a+bx+cx^2$ dont l'unique extremum est atteint lorsque la pente est nulle).
 $$
 \begin{align*}
@@ -188,7 +198,8 @@ Pour un polynôme de degré $N-1$, les fonctions de bases mentionnées ci-dessus
 La fonction ci-dessous résout le système linéaire correspondant à la matrice de Gram pour un polynôme de degré fixé. Elle retourne les coefficients du polynôme résultat.
 ```{r}
 polyreg2 <- function(data, degre) {
-  A = outer(data$X, 0:degre, "^")
+  xs <- c(data$X)
+  A = outer(xs, 0:degre, "^")
   gram = t(A) %*% A
   solve(gram, as.vector(t(A) %*% data$Y))
 }
@@ -226,4 +237,103 @@ $$
 avec \quad 0 \leq \alpha
 \end{align*}
 $$
-Il s'agit à nouveau d'un problème de minimisation quadratique en $\vec{x}$ dont le minimum peut être découvert par annulation de la dérivée.
+Il s'agit à nouveau d'un problème de minimisation quadratique en $\vec{x}$ dont le minimum peut être découvert par annulation de la dérivée.
+$$
+\begin{align*}
+ & \vec{0} = 2\vec{A}^T\vec{A}\vec{x} - 2\vec{A}^T\vec{b} + 2\alpha\vec{x} \\
+=& \\
+ & \left( \vec{A}^T\vec{A} + \alpha \vec{I}_{n\times n} \right) \vec{x} = \vec{A}^T \vec{b}
+\end{align*}
+$$
+En pratique, il s'agit donc d'ajouter une petite valeur positive $\alpha$ aux éléments de la diagonale de la matrice de Gram. Cette approche est parfois appelée "régularisation de Tikhonov" ou encore "régression Ridge".
+```{r}
+ridge <- function(alpha, data, degre) {
+  xs <- c(data$X)
+  A <- outer(xs, 0:degre, "^")
+  gram <- t(A) %*% A
+  diag(gram) <- diag(gram) + alpha
+  solve(gram, as.vector(t(A) %*% data$Y))
+}
+```
+Sur notre exemple synthétique, on affiche la fonction génératrice, le jeu de donnée et le polynôme de degré au plus sept découvert par régression ridge avec une valeur de $\alpha$ égale soit à $0$, soit à $10^{-4}$, soit à $1$.
+```{r}
+par(mfrow=c(1,3))
+coef <- ridge(0, data, 7)
+plt(data,f,main=expression(paste(plain("Degré = "), 7, plain(", "), alpha, plain(" = 0"))))
+pltpoly(coef)
+coef <- ridge(1E-4, data, 7)
+plt(data,f,main=expression(paste(plain("Degré = "), 7, plain(", "), alpha, plain(" = 1E-4"))))
+pltpoly(coef)
+coef <- ridge(1, data, 7)
+plt(data,f,main=expression(paste(plain("Degré = "), 7, plain(", "), alpha, plain(" = 1"))))
+pltpoly(coef)
+```
+Plus le coefficient de régularisation $\alpha$ est faible, moins il contraint les paramètres (c'est-à-dire, coefficients du polynômes) à conserver de petites valeurs et plus le polynôme découvert peut être complexe, au risque de provoquer du surapprentissage et donc de limiter la capacité du modèle à bien prédire les étiquettes de nouvelles observations. A l'inverse, plus le coefficient de régularisation est élevé, plus le modèle découvert sera simple, au risque de sousapprendre en ne modélisant pas convenablement les variations intrinsèques au processus qui a généré les observations.
+
+Pour mieux se représenter l'effet du coefficient de régularisation ridge, on affiche les valeurs des coefficients du polynôme découvert pour différentes valeurs de $\alpha$. Plus $\alpha$ augmente, plus les coefficients du polynôme diminuent et tendent vers $0$.
+```{r}
+alphas <- c(1E-5, 1E-4, 1E-3, 1E-2, 1E-1, 1)
+lcoef <- sapply(alphas, ridge, data, 7)
+matplot(alphas, t(lcoef), type=c("b"), pch=1, col=1:8, log="x")
+legend("topright", legend = 0:7, col=1:8, pch=1)
+```
+# Validation croisée
+Comment choisir la valeur du coefficient de régularisation $\alpha$ pour une régression ridge ?
+
+Une possibilité est de diviser le jeu de données en deux parties, l'une utilisée pour apprendre le modèle prédictif, l'autre utilisée pour valider la qualité des prédictions sur des données qui n'ont pas été vues pendant la phase d'apprentissage. On parle de jeu d'entraînement et de jeu de validation ou de test. Cette méthode est appelée "validation croisée".
+
+L'idée est de tester plusieurs valeurs de l'hyperparamètre $\alpha$ sur le jeu d'entraînement et de conserver celle qui donne les meilleurs résultats sur le jeu de test.
+
+On écrit une fonction qui permet de diviser le jeu de données en un jeu d'entraînement et un jeu de validation. On lui passe en paramètre le jeu de données initial et la proportion des données conservées pour l'entraînement.
+```{r}
+xvalpart <- function(data,p) {
+  n <- nrow(data$X)
+  nentr <- round(p*n)
+  entridx <- sample(1:n, nentr, replace=FALSE)
+  list(entr  = list(X = data$X[entridx,,drop=FALSE],  Y = data$Y[entridx]),
+       valid = list(X = data$X[-entridx,,drop=FALSE], Y = data$Y[-entridx]))
+}
+```
+
+On écrit ensuite une fonction qui apprend plusieurs modèles de régression ridge sur le jeu d'entraînement pour des valeurs différentes de l'hyperparamètre $\alpha$ et qui retourne les coefficients du meilleur modèle avec la moyenne de la valeur absolue des erreurs commises par ce modèle sur le jeu de validation.
+```{r}
+xvalridge <- function(alphas, data, degre, p) {
+  tmp <- xvalpart(data,p)
+  lcoef <- sapply(alphas, ridge, tmp$entr, degre)
+  pred <- sapply(split(lcoef,col(lcoef)), polyeval, tmp$valid$X)
+  meanabserrs <- colMeans(abs(pred - tmp$valid$Y))
+  minmeanabserr <- min(meanabserrs)
+  minidx <- which(meanabserrs == minmeanabserr)
+  list(meanabserr = minmeanabserr, coef = lcoef[,minidx], alpha = alphas[minidx])
+}
+```
+On génère un nouveau jeu de données composé de 100 observations et on calcule par validation croisée un polynôme de degré au plus égal à 5 qui modélise au mieux ces données. 
+```{r}
+data = gendat(100,0.2)
+xvalres <- xvalridge(alphas, data, 5, 0.8)
+plt(data,f)
+pltpoly(xvalres$coef)
+```
+
+# Décomposition en valeurs singulières (SVD)
+
+On présente une méthode générique de compression d'un tableau rectangulaire de données. Elle possède de nombreuses applications. En particulier, elle permet de calculer en un seul entraînement les régressions ridge pour toutes les valeurs possibles de l'hyperparamètre $\alpha$.
+
+Soit un tableau de données $\vec{X}$ composé de $N$ observations en lignes, chacune décrite par $P$ variables en colonnes. $x_{ij}$ est la valeur de la variable $j$ pour l'observation $i$. Les vecteurs lignes forment $N$ points de l'espace $\mathbb{R}^P$. Les vecteurs colonnes forment $P$ points de l'espace $\mathbb{R}^N$. On cherche à projeter le nuage des points lignes sur un sous-espace $\mathcal{H} \subset \mathbb{R}^P$, avec le moins de déformations possible. 
+
+Pour commencer, on considère le meilleur sous-espace $\mathcal{H}$ à une dimension, c'est-à-dire une droite définie par son vecteur directeur unitaire $\vec{v}$ ("unitaire" signifie que sa norme est égale à $1$, soit $\vec{v}^T\vec{v}=1$). Soit $M_i$ un des $N$ points de $\mathbb{R}^P$. A ce point correspond le vecteur $\vec{OM_i}$ aussi noté $\vec{x_i}$ car ses coordonnées sont données par la i-ème ligne de $\vec{X}$. Soit $H_i$ la projection de $M_i$ sur la droite $\mathcal{H}$.
+
+![svd1](images/svd1.jpeg)
+La longueur $OH_i$ est le produit scalaire de $\vec{x_i} et de $\vec{v}$.
+$$
+OH_i = \vec{x_i}\vec{v} = \sum_{j=1}^{P} x_{ij} v_j
+$$
+On obtient un vecteur des $N$ projections $OH_i$ par le produit de la matrice $\vec{X}$ et du vecteur $\vec{v}$ : $\vec{X}\vec{v}$. Si on choisit pour $\mathcal{H}$ le sous-espace qui minimise la somme des carrés des écarts, on doit minimiser $\sum_iM_iH_i^2$.
+
+Par le théorème de Pythagore, on a $\sum_i M_iH_i^2 = \sum_i OM_i^2 - \sum_i OH_i^2$.
+
+Puisque $\sum_i OM_i^2$ est indépendant de $\vec{v}$, on doit maximiser
+$$
+\sum OH_i^2 = (\vec{X}\vec{v})^T(\vec{X}\vec{v}) = \vec{v}^T \vec{X}^T \vec{X} \vec{v}
+$$
+...sous la contrainte $\vec{v}^T\vec{v}=1$.