Progrès...

This commit is contained in:
Pierre-Edouard Portier 2021-06-09 08:57:56 +02:00
parent 2a14a789c8
commit 4de134c21f
2 changed files with 250 additions and 73 deletions

View File

@ -1,22 +1,25 @@
---
title: "Régression"
title: "Régression et classification"
output:
pdf_document:
number_section: yes
html_notebook:
number_sections: yes
pdf_document: default
toc: yes
classoption: fleqn
---
\renewcommand{\vec}[1]{\mathbf{#1}}
# Introduction
# Concepts de base
## Premiers concepts
$\vec{x}^{(1)} \dots \vec{x}^{(P)}$ sont des vecteurs de $\mathbb{R}^N$ associés aux valeurs, aussi appelées étiquettes, $y^{(1)} \dots y^{(P)}$ de $\mathbb{R}$. On cherche une fonction $f(\vec{x}) : \mathbb{R}^N \rightarrow \mathbb{R}$ qui modélise la relation entre les observations $\vec{x}$ et les étiquettes $y$.
$\mathbf{x}^{(1)} \dots \mathbf{x}^{(P)}$ sont des vecteurs de $\mathbb{R}^N$ associés aux valeurs, aussi appelées étiquettes, $y^{(1)} \dots y^{(P)}$ de $\mathbb{R}$. Nous cherchons une fonction $f(\mathbf{x}) : \mathbb{R}^N \rightarrow \mathbb{R}$ qui modélise la relation entre les observations $\mathbf{x}$ et les étiquettes $y$.
La fonction $f$ peut avoir une forme paramétrique, comme par exemple :
$$
f(\vec{x}) = a_0 + a_1x_1 + a_2x_2 + \dots + a_Nx_N
f(\mathbf{x}) = a_0 + a_1x_1 + a_2x_2 + \dots + a_Nx_N
$$
Si $P=N+1$, les paramètres $a_0, a_1, \dots a_N$ sont solutions d'un système linéaire :
Si $P=N+1$, les paramètres $a_0, a_1, \dots a_N$ sont solution d'un système linéaire :
$$
\begin{cases}
y^{(1)} &= a_0 + a_1 x_1^{(1)} + a_2 x_2^{(1)} + \dots + a_N x_N^{(1)} \\
@ -41,15 +44,15 @@ a_0 \\ a_1 \\ \dots \\ a_N
y^{(1)} \\ y^{(2)} \\ \dots \\ y^{(P)}
\end{array} \right)
$$
Chaque ligne $i$ de la matrice du terme de gauche de l'égalité ci-dessus est le vecteur ligne $\vec{x}^{(i)T}$ avec l'addition d'un premier terme constant qui correspond au paramètre $a_0$. En nommant cette matrice $\vec{X}^T$, le système linéaire ci-dessus s'écrit :
Chaque ligne $i$ de la matrice du terme de gauche de l'égalité ci-dessus est le vecteur ligne $\mathbf{x}^{(i)T}$ avec l'addition d'un premier terme constant qui correspond au paramètre $a_0$. En nommant cette matrice $\mathbf{X}^T$, le système linéaire ci-dessus s'écrit :
$$
\vec{X}^T \vec{a} = \vec{y}
\mathbf{X}^T \mathbf{a} = \mathbf{y}
$$
On considère le cas particulier où $x$ est un scalaire et $f$ est un polynome de degré $N$ :
Soit le cas particulier où $x$ est un scalaire et $f$ est un polynome de degré $N$ :
$$
f(x) = a_0 + a_1x + a_2x^2 + \dots + a_Nx^N
$$
Avec $P=N+1$ observations et les étiquettes associées $\left( x^{(k)}, y^{(k)} \right)$, les coefficients de ce polynome sont solution d'un système linéaire :
Avec $P=N+1$ observations et les étiquettes associées $\left( x^{(k)}, y^{(k)} \right)$, les coefficients de ce polynôme sont solution d'un système linéaire :
$$
\left( \begin{array}{ccccc}
@ -69,13 +72,13 @@ $$
La matrice du terme de gauche de l'égalité ci-dessus est traditionnellement appelée "matrice de Vandermonde".
# Exemple sur un jeu de données synthétique
## Interpolation polynomiale sur un jeu de données synthétique
On considère un exemple de fonction non-linéaire.
Soit un exemple de fonction non-linéaire :
```{r}
f <- function(x) {exp(x) * cos(2*pi*sin(pi*x))}
```
On utilise cette fonction pour générer un jeu de données avec l'ajout d'un bruit gaussien.
Cette fonction est utilisée pour générer un jeu de données avec l'ajout d'un bruit gaussien.
```{r}
gendat <- function(N, sd) {
# N: nombre d'observations à générer
@ -89,7 +92,7 @@ gendat <- function(N, sd) {
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.
Voici 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, ...) {
xs = seq(0,1,length.out=100)
@ -97,7 +100,7 @@ plt <- function(data, f, ...) {
points(data$X, data$Y)
}
```
On affiche un exemple de jeu de données.
Nous affichons un exemple de jeu de données.
```{r}
set.seed(1123)
data = gendat(10,0.2)
@ -127,27 +130,27 @@ pltpoly <- function(coef) {
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.
Nous affichons la fonction génératrice, le jeu de donnée et le polynôme qui passe par chaque point du jeu de données.
```{r}
coef = polyreg1(data)
plt(data,f)
pltpoly(coef)
```
Ce polynôme qui passe exactement par chaque observation a peu de chance d'offrir de bonnes capacités prédictives. Vérifier par exemple que, sur notre exemple synthétique, pour cinq points générés à partir de la fonction $f$ et avec l'ajout d'un bruit gaussien (par exemple d'écart type $0.2$), le polynôme découvert, de degré quatre, peut être très éloigné de la fonction génératrice. C'est un exemple du phénomène de sur-apprentissage. Pour limiter ce problème, on cherche à découvrir un polynôme de degré plus faible, qui ne passera donc pas exactement par toutes les observations, mais qui aura probablement une meilleure capacité à prédire les étiquettes de nouvelles observations.
Il est improbable que ce polynôme, passant exactement par chaque observation, puisse offrir de bonnes capacités prédictives. Vérifier par exemple que, sur notre exemple synthétique, pour cinq points générés à partir de la fonction $f$ et avec l'ajout d'un bruit gaussien (par exemple d'écart type $0.2$), le polynôme découvert, de degré quatre, peut être très éloigné de la fonction génératrice. C'est un exemple du phénomène de sur-apprentissage. Pour limiter ce problème, nous cherchons à découvrir un polynôme de degré plus faible. Il ne passera pas exactement par toutes les observations, mais il prédira probablement mieux les étiquettes associées à de nouvelles observations.
# Généralisation à un espace de fonctions
On considère un espace vectoriel composé de fonctions. Une base de cet espace est un ensemble de fonctions ($f_1, f_2, \dots f_N$) tel que toute fonction de l'espace s'exprime comme combinaison linéaire des fonctions de base.
Soit un espace vectoriel composé de fonctions. Une base de cet espace est un ensemble de fonctions ($f_1, f_2, \dots f_N$) tel que toute fonction de l'espace s'exprime comme combinaison linéaire des fonctions de base.
$$
f(\vec{x}) = a_1 f_1(\vec{x}) + a_2 f_2(\vec{x}) + \dots + a_N f_N(\vec{x})
f(\mathbf{x}) = a_1 f_1(\mathbf{x}) + a_2 f_2(\mathbf{x}) + \dots + a_N f_N(\mathbf{x})
$$
Pour un jeu de données $\left\{ \vec{x^{(k)}},y^{(k)}\right\}_{k=1}^{N}$ de taille $N$, les coefficients $a_i$ sont solutions d'un système linéaire.
Pour un jeu de données $\left\{ \mathbf{x^{(k)}},y^{(k)}\right\}_{k=1}^{N}$ de taille $N$, les coefficients $a_i$ sont solution d'un système linéaire.
$$
\left( \begin{array}{ccccc}
f_1(\vec{x^{(1)}}) & f_2(\vec{x^{(1)}}) & \dots & f_N(\vec{x^{(1)}}) \\
f_1(\vec{x^{(2)}}) & f_2(\vec{x^{(2)}}) & \dots & f_N(\vec{x^{(2)}}) \\
f_1(\mathbf{x^{(1)}}) & f_2(\mathbf{x^{(1)}}) & \dots & f_N(\mathbf{x^{(1)}}) \\
f_1(\mathbf{x^{(2)}}) & f_2(\mathbf{x^{(2)}}) & \dots & f_N(\mathbf{x^{(2)}}) \\
\dots & \dots & \dots & \dots \\
f_1(\vec{x^{(N)}}) & f_2(\vec{x^{(N)}}) & \dots & f_N(\vec{x^{(N)}})
f_1(\mathbf{x^{(N)}}) & f_2(\mathbf{x^{(N)}}) & \dots & f_N(\mathbf{x^{(N)}})
\end{array} \right)
\left( \begin{array}{c}
a_1 \\ a_2 \\ \dots \\ a_N
@ -157,44 +160,44 @@ a_1 \\ a_2 \\ \dots \\ a_N
y^{(1)} \\ y^{(2)} \\ \dots \\ y^{(N)}
\end{array} \right)
$$
On note ce système linéaire $\vec{Ax} = \vec{b}$.
Nous notons ce système linéaire $\mathbf{Ax} = \mathbf{b}$.
# Approche des moindres carrés
Le système linéaire $\vec{Ax} = \vec{b}$ avec $\vec{A} \in \mathbb{R}^{M \times N}$ n'a pas de solution quand le nombre d'observations dépasse le nombre de fonctions de base (c'est-à-dire, $M>N$). Une approche est alors de chercher une approximation $\vec{Ax} \approx \vec{b}$ qui minimise la somme des carrés des erreurs : $\|\vec{A}\vec{x}-\vec{b}\|^2_2$.
## Expression matricielle
Le système linéaire $\mathbf{Ax} = \mathbf{b}$ avec $\mathbf{A} \in \mathbb{R}^{M \times N}$ n'a pas de solution quand le nombre d'observations dépasse le nombre de fonctions de base (c'est-à-dire, $M>N$). Une approche possible est alors de chercher une approximation $\mathbf{Ax} \approx \mathbf{b}$ qui minimise la somme des carrés des erreurs : $\|\mathbf{Ax}-\mathbf{b}\|^2_2$.
$$
\begin{align*}
& \|\vec{A}\vec{x}-\vec{b}\|^2_2 \\
= \{ & \|\vec{x}\|_2 = \sqrt{\vec{x}\cdot\vec{x}} \} \\
& \left(\vec{A}\vec{x}-\vec{b}\right) \cdot \left(\vec{A}\vec{x}-\vec{b}\right) \\
& \|\mathbf{A}\mathbf{x}-\mathbf{b}\|^2_2 \\
= \{ & \|\mathbf{x}\|_2 = \sqrt{\mathbf{x}\cdot\mathbf{x}} \} \\
& \left(\mathbf{A}\mathbf{x}-\mathbf{b}\right) \cdot \left(\mathbf{A}\mathbf{x}-\mathbf{b}\right) \\
= \{ & \text{Par définition du produit scalaire euclidien} \} \\
& \left(\vec{A}\vec{x}-\vec{b}\right)^T \left(\vec{A}\vec{x}-\vec{b}\right) \\
& \left(\mathbf{A}\mathbf{x}-\mathbf{b}\right)^T \left(\mathbf{A}\mathbf{x}-\mathbf{b}\right) \\
= \{ & \text{propriété de la transposition} \} \\
& \left(\vec{x}^T\vec{A}^T - \vec{b}^T \right) \left(\vec{A}\vec{x}-\vec{b}\right) \\
& \left(\mathbf{x}^T\mathbf{A}^T - \mathbf{b}^T \right) \left(\mathbf{A}\mathbf{x}-\mathbf{b}\right) \\
= \{ & \text{multiplication} \} \\
& \vec{x}^T\vec{A}^T\vec{A}\vec{x} - \vec{x}^T\vec{A}^T\vec{b} - \vec{b}^T\vec{A}\vec{x} + \vec{b}^T\vec{b} \\
= \{ & \vec{b}^T\vec{A}\vec{x} \text{ étant une valeur scalaire, } \vec{b}^T\vec{A}\vec{x} = \left(\vec{b}^T\vec{A}\vec{x}\right)^T = \vec{x}^T\vec{A}^T\vec{b} \} \\
& \vec{x}^T\vec{A}^T\vec{A}\vec{x} - 2\vec{x}^T\vec{A}^T\vec{b} + \vec{b}^T\vec{b}
& \mathbf{x}^T\mathbf{A}^T\mathbf{A}\mathbf{x} - \mathbf{x}^T\mathbf{A}^T\mathbf{b} - \mathbf{b}^T\mathbf{A}\mathbf{x} + \mathbf{b}^T\mathbf{b} \\
= \{ & \mathbf{b}^T\mathbf{A}\mathbf{x} \text{ étant une valeur scalaire, } \mathbf{b}^T\mathbf{A}\mathbf{x} = \left(\mathbf{b}^T\mathbf{A}\mathbf{x}\right)^T = \mathbf{x}^T\mathbf{A}^T\mathbf{b} \} \\
& \mathbf{x}^T\mathbf{A}^T\mathbf{A}\mathbf{x} - 2\mathbf{x}^T\mathbf{A}^T\mathbf{b} + \mathbf{b}^T\mathbf{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).
$$
Cette dernière expression quadratique en $\mathbf{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*}
& \vec{0} = 2\vec{A}^T\vec{A}\vec{x} - 2\vec{A}^T\vec{b} \\
& \mathbf{0} = 2\mathbf{A}^T\mathbf{A}\mathbf{x} - 2\mathbf{A}^T\mathbf{b} \\
=& \\
& \vec{A}^T\vec{A}\vec{x} = \vec{A}^T\vec{b}
& \mathbf{A}^T\mathbf{A}\mathbf{x} = \mathbf{A}^T\mathbf{b}
\end{align*}
Ainsi, quand $M>N$, la solution approximée $\mathbf{x}$, telle que $\mathbf{Ax} \approx \mathbf{b}$ par minimisation de la somme des carrés des erreurs, est la solution du système linéaire suivant où $\mathbf{A}^T\mathbf{A}$ est appelée la matrice de Gram.
$$
Ainsi, quand $M>N$, la solution approximée $\vec{x}$, telle que $\vec{Ax} \approx \vec{b}$ en minimisant la somme des carrés des erreurs, est la solution du système linéaire suivant où $\vec{A}^T\vec{A}$ est appelée la matrice de Gram.
$$
\vec{A}^T\vec{A} \vec{x} = \vec{A}^T\vec{b}
\mathbf{A}^T\mathbf{A} \mathbf{x} = \mathbf{A}^T\mathbf{b}
$$
## Approche des moindres carrés appliquée à la régression polynomiale
Pour un polynôme de degré $N-1$, les fonctions de bases mentionnées ci-dessus sont : $f_1(x)=1$, $f_2(x)=x$, $f_3(x)=x^2$,..., $f_N(x)=x^{N-1}$. Elles permettent de définir la matrice $\vec{A}$ et la matrice de Gram $\vec{A}^T\vec{A}$.
Pour un polynôme de degré $N-1$, les fonctions de bases mentionnées ci-dessus sont : $f_1(x)=1$, $f_2(x)=x$, $f_3(x)=x^2$,..., $f_N(x)=x^{N-1}$. Elles permettent de définir la matrice des données $\mathbf{A}$ et la matrice de Gram $\mathbf{A}^T\mathbf{A}$.
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) {
@ -204,7 +207,7 @@ polyreg2 <- function(data, degre) {
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 "meilleur" polynôme de degré 3.
Sur notre exemple synthétique, nous affichons la fonction génératrice, le jeu de donnée et le "meilleur" polynôme de degré 3.
```{r}
coef = polyreg2(data,3)
plt(data,f)
@ -214,7 +217,7 @@ Ce polynôme de degré trois modélise mieux la fonction génératrice inconnue
# Régularisation de Tikhonov
Avec moins d'observations que de fonctions de base ($M<N$), le système $\vec{A}\vec{x}=\vec{b}$ ne possède pas de solution unique. Même quand $M \geq N$, le système linéaire peut posséder une solution approchée préférable à la solution optimale. C'est en particulier vrai quand plusieurs observations sont très proches. Par exemple, soit le système linéaire suivant :
Avec moins d'observations que de fonctions de base ($M<N$), le système $\mathbf{A}\mathbf{x}=\mathbf{b}$ ne possède pas de solution unique. Même quand $M \geq N$, le système linéaire peut posséder une solution approchée préférable à la solution optimale. C'est en particulier vrai quand plusieurs observations sont très proches à un coefficient multiplicatif près (on parle de colinéarités). Par exemple, soit le système linéaire suivant :
$$
\left( \begin{array}{cc}
1 & 1 \\
@ -228,24 +231,24 @@ x_1 \\ x_2
1 \\ 0.99
\end{array} \right)
$$
Sa solution est $\vec{x}^T = (1001,-1000)$. Cependant, la solution approchée $\vec{x}^T = (0.5,0.5)$ est nettement préférable. En effet, la solution optimale a peu de chance de bien s'adapter à de nouvelles observations (par exemple, l'observation $(1,2)$ serait projetée sur l'étiquette -999).
Sa solution est $\mathbf{x}^T = (1001,-1000)$. Cependant, la solution approchée $\mathbf{x}^T = (0.5,0.5)$ semble préférable. En effet, la solution optimale a peu de chance de bien s'adapter à de nouvelles observations (par exemple, l'observation $(1,2)$ serait projetée sur l'étiquette -999).
Ainsi, lorsqu'il faut choisir entre plusieurs solutions, il peut être efficace d'exprimer une préférence envers celles dont les coefficients (ou paramètres) ont de faibles valeurs. Cela consiste par exemple à minimiser $|x_1|+|x_2|+\dots$ (aussi noté $\|\mathbf{x}\|_1$, la "norme 1") ou encore $x_1^2+x_2^2+\dots$ (aussi noté $\|\mathbf{x}\|_2^2$, le carré de la "norme 2"). Dans ce dernier cas, il s'agit de résoudre un nouveau problème de minimisation :
Ainsi, lorsqu'il faut choisir entre plusieurs solutions, on préférera celles dont les coefficients (ou paramètres) ont de faibles valeurs. Par exemple, on pourra vouloir minimiser $|x_1|+|x_2|+\dots$ (aussi noté $\|\vec{x}\|_1$, la "norme 1") ou encore $x_1^2+x_2^2+\dots$ (aussi noté $\|\vec{x}\|_2^2$, le carré de la "norme 2"). Dans ce dernier cas, on peut chercher à résoudre un nouveau problème de minimisation :
$$
\begin{align*}
\min_{\vec{x}} \|\vec{A}\vec{x}-\vec{b}\|^2_2 + \alpha \|\vec{x}|^2_2 \\
\min_{\mathbf{x}} \|\mathbf{A}\mathbf{x}-\mathbf{b}\|^2_2 + \alpha \|\mathbf{x}|^2_2 \\
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.
$$
C'est encore un problème de minimisation quadratique en $\mathbf{x}$ dont le minimum se découvre 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} \\
& \mathbf{0} = 2\mathbf{A}^T\mathbf{A}\mathbf{x} - 2\mathbf{A}^T\mathbf{b} + 2\alpha\mathbf{x} \\
=& \\
& \left( \vec{A}^T\vec{A} + \alpha \vec{I}_{n\times n} \right) \vec{x} = \vec{A}^T \vec{b}
& \left( \mathbf{A}^T\mathbf{A} + \alpha \mathbf{I}_{n\times n} \right) \mathbf{x} = \mathbf{A}^T \mathbf{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".
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 porte plusieurs noms dont "régularisation de Tikhonov" ou "régression Ridge".
```{r}
ridge <- function(alpha, data, degre) {
xs <- c(data$X)
@ -255,7 +258,7 @@ ridge <- function(alpha, data, degre) {
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$.
Sur notre exemple synthétique, nous affichons 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)
@ -268,23 +271,25 @@ 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.
Plus le coefficient de régularisation $\alpha$ est faible, moins il contraint les paramètres (c'est-à-dire, les coefficients du polynômes) à conserver de petites valeurs et plus le polynôme découvert peut être complexe, au risque de provoquer du sur-apprentissage 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 sous-apprendre en ne modélisant pas convenablement les variations propres 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$.
Pour mieux visualiser l'effet du coefficient de régularisation ridge, nous affichons 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.
Nous écrivons une fonction qui permet de diviser le jeu de données en un jeu d'entraînement et un jeu de validation. Nous lui passons 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)
@ -295,7 +300,7 @@ xvalpart <- function(data,p) {
}
```
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.
Nous écrivons 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)
@ -307,7 +312,7 @@ xvalridge <- function(alphas, data, degre, p) {
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.
Nous générons un nouveau jeu de données composé de 100 observations et nous calculons 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)
@ -317,23 +322,195 @@ 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$.
## Formulation du SVD comme un problème d'optimisation
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.
Nous présentons 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$.
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}$.
Soit un tableau de données $\mathbf{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$. Nous cherchons à projeter le nuage des points lignes sur un sous-espace $\mathcal{H} \subset \mathbb{R}^P$, tout en minimisant les déformations.
Pour commencer, nous considèrons le meilleur sous-espace $\mathcal{H}$ à une dimension, c'est-à-dire une droite définie par son vecteur directeur unitaire $\mathbf{v}$ ("unitaire" signifie que sa norme est égale à $1$, soit $\mathbf{v}^T\mathbf{v}=1$). Soit $M_i$ un des $N$ points de $\mathbb{R}^P$. A ce point correspond le vecteur $\mathbf{OM_i}$ aussi noté $\mathbf{x_i}$ car ses coordonnées sont données par la i-ème ligne de $\mathbf{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}$.
La longueur $OH_i$ est le produit scalaire de $\mathbf{x_i}$ et de $\mathbf{v}$.
$$
OH_i = \vec{x_i}\vec{v} = \sum_{j=1}^{P} x_{ij} v_j
OH_i = \mathbf{x_i}\mathbf{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$.
Nous obtenons un vecteur des $N$ projections $OH_i$ par le produit de la matrice $\mathbf{X}$ et du vecteur $\mathbf{v}$ : $\mathbf{X}\mathbf{v}$. En choisissant pour $\mathcal{H}$ le sous-espace qui minimise la somme des carrés des écarts, il faut 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$.
Le théorème de Pythagore donne $\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
Puisque $\sum_i OM_i^2$ est indépendant de $\mathbf{v}$, nous devons maximiser
$$
\sum OH_i^2 = (\vec{X}\vec{v})^T(\vec{X}\vec{v}) = \vec{v}^T \vec{X}^T \vec{X} \vec{v}
\sum OH_i^2 = (\mathbf{X}\mathbf{v})^T(\mathbf{X}\mathbf{v}) = \mathbf{v}^T \mathbf{X}^T \mathbf{X} \mathbf{v}
$$
...sous la contrainte $\vec{v}^T\vec{v}=1$.
...sous la contrainte $\mathbf{v}^T\mathbf{v}=1$.
## Matrices définies non négatives
Avant de résoudre ce problème d'optimisation, nous rappelons quelques propriétés des matrices symétriques définies non négatives. Une matrice définie non négative (souvent noté PSD pour "positive semi-definite") $\mathbf{M}$ représente une métrique car elle définit un produit scalaire et donc une géométrie, c'est-à-dire qu'elle donne un sens aux notions de distance et d'angle.
Le produit scalaire défini par $\mathbf{M}$ entre deux vecteurs $\mathbf{x}$ et $\mathbf{y}$ pourra se noter $<\mathbf{x},\mathbf{y}>_{\mathbf{M}} = \mathbf{x}^T\mathbf{M}\mathbf{y}$. Pour que ce produit scalaire soit valide, $\mathbf{M}$ doit respecter la contrainte $\mathbf{x}^T\mathbf{M}\mathbf{y} \geq 0$ pour tout $\mathbf{x}$ et $\mathbf{y}$.
Une fois donné un produit scalaire, les notions liées de norme et d'angle suivent~: $\|\mathbf{x}\|_{\mathbf{M}} = \sqrt{<\mathbf{x},\mathbf{x}>_{\mathbf{M}}}$ et $cos^2(\mathbf{x},\mathbf{y}) = <\mathbf{x},\mathbf{y}>_{\mathbf{M}}/\|\mathbf{x}\|_{\mathbf{M}}\|\mathbf{y}\|_{\mathbf{M}}$.
La géométrie euclidienne "classique" (associée à la base canonique) correspond à l'emploi de la matrice identité pour métrique, $\mathbf{M}=\mathbf{I}$.
Aussi, les valeurs propres d'une matrice PSD sont toutes positives~:
\begin{align*}
& \mathbf{x}^T\mathbf{M}\mathbf{x} \geq 0 \\
= \{ & \lambda \text{ est une valeur propre de } \mathbf{M} \} \\
& \mathbf{x}^T\lambda\mathbf{x} \geq 0 \\
= \{ &\text{utilisation de la métrique identité} \} \\
& \lambda \|x\|^2_{\mathbf{I}} \geq 0 \\
\Rightarrow \{ &\|x\|^2_{\mathbf{I}} \geq 0 \} \\
& \lambda \geq 0
\end{align*}
Soit $\mathbf{V}$ la matrice des vecteurs propres de $\mathbf{M}$ et $\mathbf{L}$ la matrice diagonale de ses valeurs propres positives.
\begin{align*}
& \mathbf{M}\mathbf{V} = \mathbf{V}\mathbf{L} \\
= \phantom{\{} & \\
& \mathbf{M} = \mathbf{V}\mathbf{L}\mathbf{V}^{-1} \\
= \{ & \mathbf{S} = \sqrt{\mathbf{L}}\} \\
& \mathbf{M} = \mathbf{V}\mathbf{S}\mathbf{S}\mathbf{V}^{-1} \\
= \{ & \quad \text{Les vecteurs propres distincts sont orthogonaux deux à deux : } \mathbf{V}^{-1} = \mathbf{V}^T. \\
\phantom{=}\phantom{ }\phantom{\{} & \quad \mathbf{S}=\mathbf{S}^T. \text{On pose } \mathbf{T}=\mathbf{V}\mathbf{S}. \} \\
& \mathbf{M} = \mathbf{T}\mathbf{T}^T \\
\end{align*}
Ainsi, toute matrice définie non négative $\mathbf{M}$ peut être décomposée sous la forme $\mathbf{M}=\mathbf{T}\mathbf{T}^T$ où $\mathbf{T}$ est une transformation linéaire composée d'un changement d'échelle des axes du repère ($\mathbf{S}$) et d'une rotation ($\mathbf{V}$).
Nous remarquons par exemple que l'équation d'un cercle selon la métrique $\mathbf{M}$ correspond à celle d'une ellipse selon la métrique identité $\mathbf{I}$ :
\begin{align*}
& \|\mathbf{x}\|_{\mathbf{M}}^2 = c \\
= \{ &\text{C'est l'équation d'un cercle, c est une constante.}\} \\
& < \mathbf{x}, \mathbf{x}>_{\mathbf{M}} = c \\
= \{ &\text{Par définition du produit scalaire.} \} \\
& \mathbf{x}^T \mathbf{M} \mathbf{x} = c \\
= \{ &\mathbf{M} = \mathbf{T}^T\mathbf{T}\} \\
& \mathbf{x}^T \mathbf{T}^T \mathbf{T} \mathbf{x} = c \\
= \{ &\text{Par définition du produit scalaire, on retrouve le cercle déformé en une ellipse par } \mathbf{T} \} \\
& \|\mathbf{T}\mathbf{x}\|_{\mathbf{I}}^2 = c \\
\end{align*}
Ainsi, une matrice symétrique définie non négative projette le cercle unité sur une ellipse dont les axes orthonormaux sont les vecteurs propres et les longueurs des axes sont les valeurs propres.
Enfin, l'angle formé entre $\mathbf{x}$ et $\mathbf{Mx}$ est inférieur ou égal à $\pi/2$ radians :
\begin{align*}
& \mathbf{x}^T \mathbf{M} \mathbf{x} \geq 0 \\
= \{&\text{Géométrie du produit scalaire. On note } \theta \text{ l'angle entre les deux vecteurs.}\} \\
& \|\mathbf{x}\|_\mathbf{I} \; \|\mathbf{Mx}\|_\mathbf{I} \; cos(\theta) \geq 0 \\
\Rightarrow \phantom{\{}&\\
& |\theta| \leq \pi/2
\end{align*}
Ainsi, le vecteur résultat $\mathbf{Mx}$ est contraint de rester du même côté que $\mathbf{x}$ de l'hyperplan perpendiculaire à $\mathbf{x}$ qui sépare l'espace en deux. Nous comprenons que le concept de matrice symétrique définie non négative est une généralisation pour un espace multidimensionnel du concept de nombre positif sur la droite réelle.
## Dérivation de la décomposition en valeur singulière
Après cette parenthèse sur la notion de métrique associée à une matrice définie non négative, nous reprenons la dérivation de la décomposition en valeurs singulières.
Nous notons $\mathbf{v_1}$ le vecteur directeur du meilleur sous-espace de dimension 1 pour le nuage des $N$ observations de $\mathcal{R}^P$. Le meilleur sous-espace de dimension $2$ contient $\mathbf{v_1}$ et il faut le compléter avec un second vecteur unitaire $\mathbf{v_2}$, orthogonal à $\mathbf{v_1}$ et qui maximise $\mathbf{v_2}^T \mathbf{X}^T \mathbf{X} \mathbf{v_2}$.
Nous montrons maintenant que $\mathbf{v_1}$ est le vecteur propre de $\mathbf{X}^T\mathbf{X}$ associé à la plus grande valeur propre $\lambda_1$.
Le résultat est sans difficulté plus général. Nous le prouvons pour toute matrice symétrique $\mathbf{A}$ (non seulement $\mathbf{X}^T\mathbf{X}$) et toute métrique de $\mathcal{R}^P$ représentée par une matrice symétrique définie non négative $\mathbf{M}$ (non seulement la martice identité $\mathbf{I}$).
Nous rappelons que $\mathbf{v}^T\mathbf{A}\mathbf{v} = \Sigma a_{i,j}v_iv_j$. Ainsi, comme $\mathbf{A}$ et $\mathbf{M}$ sont symétriques (i.e., $a_{i,j}=a_{j,i}$ et $m_{i,j}=m_{j,i}$), les dérivées partielles des formes quadratiques ont les formes suivantes :
$$
\frac{\partial \mathbf{v}^T\mathbf{A}\mathbf{v}}{\partial \mathbf{v}}=2\mathbf{A}\mathbf{v} \quad \text{et} \quad \frac{\partial \mathbf{v}^T\mathbf{M}\mathbf{v}}{\partial \mathbf{v}}=2\mathbf{M}\mathbf{v}
$$
Pour trouver le maximum de $\mathbf{v}^T\mathbf{A}\mathbf{v}$ en intégrant la contrainte unitaire sur $\mathbf{v}$ (i.e., $\mathbf{v}^T\mathbf{M}\mathbf{v}=1$), nous annulons les dérivées du langrangien $\mathcal{L}$ :
$$
\mathcal{L} = \mathbf{v}^T\mathbf{A}\mathbf{v} - \lambda (\mathbf{v}^T\mathbf{M}\mathbf{v} - 1)
$$
\begin{align*}
& \frac{\partial \mathcal{L}}{\partial \mathbf{v}} = 0 \\
= \{& \text{voir calcul des dérivées ci-dessus} \} \\
& 2\mathbf{A}\mathbf{v} - 2\lambda\mathbf{M}\mathbf{v} = 0 \\
= \phantom{\{} & \\
& \mathbf{A}\mathbf{v} = \lambda\mathbf{M}\mathbf{v} \\
= \{& \mathbf{v}^T\mathbf{M}\mathbf{v}=1 \} \\
& \lambda = \mathbf{v}^T\mathbf{A}\mathbf{v}
\end{align*}
Nous découvrons ainsi que la valeur du multiplicateur de Lagrange $\lambda$ est le maximum recherché. Par ailleurs, nous avons :
\begin{align*}
& \mathbf{A}\mathbf{v} = \lambda\mathbf{M}\mathbf{v} \\
= \{& \mathbf{M} \text{ est définie non négative et donc inversible} \} \\
& \mathbf{M}^{-1}\mathbf{A}\mathbf{v} = \lambda\mathbf{v}
\end{align*}
Donc $\mathbf{v}$ est le vecteur propre de $\mathbf{M}^{-1}\mathbf{A}$ associé à la plus grande valeur propre $\lambda$. Nous notons cette solution $\mathbf{v_1}$ et la valeur propre correspondante $\lambda_1$.
Nous cherchons ensuite $\mathbf{v_2}$, unitaire ($\mathbf{v_2}^T\mathbf{M}\mathbf{v_2}=1$), orthogonal à $\mathbf{v_1}$ ($\mathbf{v_1}^T\mathbf{M}\mathbf{v_2}=0$) et qui maximise $\mathbf{v_2}^T\mathbf{A}\mathbf{v_2}$. Pour ce faire, nous annulons les dérivées du Langrangien ci-après.
$$
\mathcal{L} = \mathbf{v_2}^T\mathbf{A}\mathbf{v_2} - \lambda_2 (\mathbf{v_2}^T\mathbf{M}\mathbf{v_2} - 1) - \mu_2\mathbf{v_2}^T\mathbf{M}\mathbf{v_1}
$$
En cours de rédaction...
# Régression avec régularisation de Tikhonov et SVD
# Validation croisée "un contre tous" (leave-one-out-cross-validation, LOOCV) et régression régularisée
# Interprétation d'un modèle de régression linéraire
## Signification des coefficients
## Exemple sur un dataset
## Intervalles de confiances et approche par bootstrap
# Estimateurs par maximum de vraissemblance
## Application aux paramètres d'une loi normale uni-dimensionnelle
Permet d'introduire les notions de biais et de variance minimale.
## Notation généralisée pour toute distribution
## Vraissemblance conditionnelle pour la régression linéaire
## Interprétation bayésienne de la régularisation commme distribution a priori des paramètres
# Régularisation et sélection de variable avec l'approche LASSO
# Régression logistique
## Vraissemblance conditionnelle pour la régression logistique
## Interprétation des coefficients d'une régression logistique
## Descente de gradient (stochastique) appliquée à la régression logistique
# Relations générales entre biais, variance, nombre d'observations et complexité du modèle
# Ensembles de modèles
## Effets sur le biais et la variance d'une approche par ensembles de modèles
## Bagging illustré avec le modèle de forêt aléatoire
## Boosting ?
# Réseaux de neurones
## Neurone
## Architecture en couches
## Rétropropagation du gradient
## Régularisation implicite par descente de gradient stochastique
## Régularisation par dropout
## Régularisation par batch normalization

BIN
ML1_regression.pdf Normal file

Binary file not shown.