Le document suivant a été réalisé avec le logiciel libre 
[https://www.r-project.org/](**R**). Les informations sur les versions précises
du logiciel et des packages ainsi que sur la configuration de la machine (en 
particulier, ce document a été réalisé sous Linux et la prise en charge des
caractères accentués est réalisée par l'encodage UTF-8 ; pour compiler 
correctement cette vignette sous Windows, il sera probablement nécessaire 
d'utiliser le menu « Fichier / Ré-ouvrir avec l'encodage ... ») au moment de la
compilation sont fournies à la fin de ce document (section « Information sur la
session **R** »).

Le but de ce document est d'illustrer, sur des cas très simples, l'utilisation
et le comportement des perceptrons à une couche cachée. Les études de cas sont
réalisées avec les packages 
[```nnet```](https://cran.r-project.org/web/packages/nnet) (apprentissage du 
réseaux de neurones) et 
[```e1071```](https://cran.r-project.org/web/packages/e1071) (calibration des
divers hyper-paramètres de la méthode). En outre, le package
[```mlbench```](https://cran.r-project.org/web/packages/mlbench) est utilisé
pour fournir un exemple de données réelles (dans la deuxième partie de ce 
document « Un exemple sur des données réelles ») et le package 
[```devtools```](https://cran.r-project.org/web/packages/devtools), non 
indispensable pour les applications, est simplement utilisé pour fournir des
informations sur le logiciel, les packages et l'environnement de compilation
utilisés.

Les packages peuvent être installés avec la commande
```{r installPkg, eval=FALSE}
install.packages(c("nnet", "e1701", "mlbench", "devtools"))
```
et ils sont chargés dans l'environnement de travail de **R** avec
```{r loadLib}
library(devtools)
library(nnet)
library(e1071)
library(mlbench)
```

Enfin, le présent document (format HTML) ainsi que sa source (format RMarkdown,
compilable en utilisant RStudio et les packages ```evaluate```, ```formatR```,
```highr```, ```markdown```, ```yaml```, ```htmltools```, ```knitr``` et 
```rmarkdown``` ; ces packages s'installent automatiquement à la première 
compilation d'un fichier au format RMarkdown sous RStudio) sont disponible sur
[mon site web](http://www.nathalievialaneix.eu/teaching/JES2016/).


# Données simulées

La première partie de ce document illustre le comportement du perceptron dans 
une tâche simple qui est l'estimation d'une fonction de R dans R à partir 
d'observations bruitées.

Les données sont générées par les commandes suivantes (dans lesquelles la 
fonction ```set.seed``` est utilisée pour initialiser une graine aléatoire et 
assurer des résultats reproductibles) :

```{r dataGeneration}
set.seed(1528)
x <- runif(15, 0, 1)
y <- sin(x*2*pi) + rnorm(15, 0, 0.3)
```

Les données générées sont représentées dans la figure ci-dessous avec la 
fonction sous-jacente à estimer, $x \rightarrow \sin(2\pi x)$ :
```{r plotData}
t <- seq(0, 1, length=1000)
yt <- sin(t*2*pi)
plot(x, y, pch = "+", lwd = 2)
lines(t, yt, col = "darkred")
```

```{r exportPlot, echo=FALSE, results='hide', fig.show='hide'}
plot(x, y, pch = "+", lwd = 2, ylim = range(c(y, yt)), cex=2, cex.axis=1.5,
     cex.lab=1.5)
lines(t, yt, col = "darkred", lwd=2)
dev.print(png, file = "../figures/simulated_data.png", width = 800, 
          height = 600, res = 150)
dev.print(pdf, file = "../figures/simulated_data.pdf", width=9.6, height=7.2)
```


## Influence de $Q$

Dans un premier temps, nous illustrons l'influence de l'hyper-paramètre $Q$ (le
nombre de neurones sur la couche cachée). L'apprentissage du perceptron à une 
couche cachée est réalisée avec la fonction ```nnet``` sur laquelle on peut 
obtenir de l'aide en utilisant la commande
```{r helpNNet, eval=FALSE}
help(nnet)
```
Dans ce qui suit, pour les observations générées précédemment et pour chaque 
valeur de $Q$ variant entre 1 et 10, 10 apprentissages sont réalisés. Le 
perceptron qui minimise l'erreur quadratique empirique sur les données 
d'apprentissage est finalement conservé (un perceptron pour chaque valeur de 
$Q$) :

```{r perceptronK}
set.seed(1513)
Qlist <- 1:10
res <- list()
iter <- 10
for (k in Qlist) {
  best_error <- 10^5
  for (ind in 1:iter) {
    cur_nnet <- nnet(x, y, size = k, trace = FALSE, linout = TRUE)
    if (mean(cur_nnet$residuals^2 < best_error))
      best_nnet <- cur_nnet
  }
  res[[k]] <- best_nnet
}
```

La fonction estimée par le perceptron est ensuite représentée en comparaison de
la fonction cible. Pour cela, on utilise la fonction ```predict``` qui permet 
d'obtenir la prédiction du perceptron pour 1000 valeurs de $x$ uniformément
réparties dans $[0,1]$. L'erreur moyenne réalisée par le perceptron sur ces 100
valeurs est conservées dans la variable ```gen_err``` : cette erreur moyenne 
sera comparée à l'erreur moyenne d'apprentissage (conservée dans la variable
```train_err```) :

```{r plotK}
par(mfrow = c(2, 5))
train_err <- unlist(lapply(res, function(anet) anet$value/length(x)))
gen_err <- vector(length = 10)
for (k in Qlist) {
  par(mar = rep(3,4))
  cur_pred <- predict(res[[k]], as.matrix(t, ncol = 1))
  plot(x, y, pch = "+", lwd = 2, main = paste0("Q = ", k),
       ylim = range(c(y, yt, cur_pred)), cex=2, cex.main=2, axes=FALSE)
  box()
  lines(t, yt, col = "darkred", lwd=2, lty=3)
  lines(t, cur_pred, col = "black", lwd=2)
  gen_err[k] <- mean((cur_pred - yt)^2)
}
```

```{r exportPlot2, echo=FALSE, results='hide', fig.show='hide'}
par(mfrow = c(2, 5))
for (k in Qlist) {
  par(mar = rep(3,4))
  cur_pred <- predict(res[[k]], as.matrix(t, ncol = 1))
  plot(x, y, pch = "+", lwd = 2, main = paste0("Q = ", k),
       ylim = range(c(y, yt, cur_pred)), cex=2, cex.main=2, axes=FALSE)
  box()
  lines(t, yt, col = "darkred", lwd=2, lty=3)
  lines(t, cur_pred, col = "black", lwd=2)
}
dev.print(png, file = "../figures/simulated_est_K.png", width = 1200, 
          height = 800, res = 150)
dev.print(pdf, file = "../figures/simulated_est_K.pdf", width = 12, height = 8)
```

Le graphique ci-dessous présente justement une comparaison de l'évolution de
l'erreur empirique d'apprentissage (```train_err```) avec l'erreur commise sur 
la grille régulière de 1000 points sur $[0,1]$ (```gen_err```) qui sera appelée
« erreur de test » et correspond à une estimation de la quantité
\[
  \mathbb{E}\left[\left(\hat{f}_Q(X) - f(X) \right)^2\right]
\]
pour $X \sim \mathcal{U}[0,1]$ et $f(x) = \sin(2\pi x)$ :

```{r plotErrK}
plot(Qlist, train_err, pch = "+", col = "darkred", type = "b", 
     ylim = c(0, max(c(train_err, gen_err))), xlab = "Q", ylab = "erreur", 
     cex.axis = 1.5, cex.lab = 1.5)
lines(Qlist, gen_err, col = "black", lty=2)
points(Qlist, gen_err, pch = "+", col = "black")
legend(0.7, 0.25, col = c("darkred", "black"), cex = 1.5, bty = "n",
       legend = c("apprentissage", "test"), lwd = 2, lty=c(1,2))
```

```{r exportPlot3, echo=FALSE, results='hide', fig.show='hide'}
plot(Qlist, train_err, pch = "+", col = "darkred", type = "b", 
     ylim = c(0, max(c(train_err, gen_err))), xlab = "Q", ylab = "erreur", 
     cex.axis = 1.5, cex.lab = 1.5)
lines(Qlist, gen_err, col = "black", lty=2)
points(Qlist, gen_err, pch = "+", col = "black")
legend(0.7, 0.25, col = c("darkred", "black"), cex = 1.5, bty = "n",
       legend = c("apprentissage", "test"), lwd = 2, lty=c(1,2))
dev.print(png, file = "../figures/erreur_K.png", width = 800, height = 600, 
          res = 150)
dev.print(pdf, file = "../figures/erreur_K.pdf", width=9.6, height=7.2)
```

On y observe que l'erreur d'apprentissage a tendance à décroître lorsque $Q$ 
augmente alors que l'erreur de test a une phase de décroissance suivi 
(grossièrement) d'une phase de croissance. Le minimum de l'erreur de test est
atteint pour $Q=$ `r which.min(gen_err)` avec une erreur égale à :

```{r minErrK}
min(gen_err)
```


## Calibration de $Q$

En l'absence d'un grand échantillon de test pour calibrer la valeur de $Q$ par
calcul d'une erreur de tests, diverses stratégies classiques (validation simple,
validation croisée, bootstrap...) peuvent être mises en œuvre pour sélectionner
une valeur pertinente pour $Q$. Ces stratégies peuvent être mise en œuvre grâce
à la fonction ```tune.nnet``` du package e1071. Les diverses options possibles
pour la calibration (type de calibration - par validation simple, croisée,
bootstrap ... - ou paramètre de la méthode utilisée comme le nombre de groupes
dans la validation croisée) sont controlées via la fonction ```tune.control```
dont le résultat est passé en argument à la fonction ```tune.nnet```.

Dans la commande ci-dessous, on illustre la callibration par validation croisée
(10 groupes) de la valeur de $Q$ pour les données précédentes. Ici, l'option
```nrepeat``` sert à fixer le nombre d'initialisations différentes (et donc de
perceptrons différents) qui sont effectuées pour produire, pour chaque groupe de
la validation croisée, une estimation différente de l'erreur :

```{r tuneK}
set.seed(1634)
res_tune <- tune.nnet(x, y, tune.control(nrepeat = iter), size = Qlist, 
                      linout = TRUE)
```

Le graphique ci-dessous illustre l'évolution de l'erreur de validation croisée
obtenue pour les différentes valeurs de $Q$ :

```{r plotTuneK}
plot(Qlist, res_tune$performances$error, main = "", xlab = "Q", ylab = "erreur",
     type = "b")
```

```{r exportPlot4, echo=FALSE, results='hide', fig.show='hide'}
plot(Qlist, res_tune$performances$error, main = "", xlab = "Q", ylab = "erreur",
     type = "b")
dev.print(png, file = "../figures/tune_K.png", width = 800, height = 600, 
          res = 150)
```

L'erreur de validation croisée est donc minimale pour $Q=$
`r which.min(res_tune$best.parameters)` et on peut obtenir les informations
additionnelles suivantes sur sa valeur et sur les caractéristiques du perceptron
obtenu à partir de ce paramètre optimal :

```{r bestNNETK}
res_tune$best.parameters
res_tune$best.performance
res_tune$best.model
```


## Régularisation

Cette partie illustre l'utilisation d'une pénalité dans l'apprentissage du
perceptron. De manière plus précise, la fonction objectif minimisée lors de 
l'apprentissage est :
\[
  \hat{\mathcal{R}}_n(w) + \lambda \|w\|^2
\]
où $\lambda$ est le paramètre de régularisation. Dans les lignes de commande 
ci-dessous, on fait varier ce paramètre selon une échelle exponentielle pour un
perceptron à une couche cachée dont le nombre de neurones est très supérieur à
l'optimal (on a choisi $Q=10$) :

```{r perceptronDecay}
set.seed(1546)
k <- 10
reg_list <- c(10^(- (1:9)), 0)
res <- list()
iter <- 10
for (ind in seq_along(reg_list)) {
  best_error <- 10^5
  for (repet in 1:iter) {
    cur_nnet <- nnet(x, y, size = k, trace = FALSE, linout = TRUE, 
                     decay = reg_list[ind])
    if (mean(cur_nnet$residuals^2 < best_error))
      best_nnet <- cur_nnet
  }
  res[[ind]] <- best_nnet
}
```

Les résultats obtenus, comme dans la partie « Calibration de $Q$ » sont 
illustrées par le graphique suivant dans lequel on constate qu'un fort paramètre
de régularisation conduit à une estimation trop grossière de la fonction cible
alors qu'un faible paramètre de régularisation produit des résultats similaires
à ceux obtenus sans régularisation, à savoir un sur-apprentissage lors de 
l'apprentissage :

```{r plotDecay}
par(mfrow = c(2,5))
train_err <- unlist(lapply(res, function(anet) anet$value/length(x)))
gen_err <- vector(length = 10)
print_reg <- c(reg_list[-10], "sans régularisation")
for (ind in seq_along(reg_list)) {
  par(mar = rep(3,4))
  cur_pred <- predict(res[[ind]], as.matrix(t, ncol = 1))
  plot(x, y, pch = "+", lwd = 2, main = print_reg[ind], 
       ylim = range(c(y, yt, cur_pred)), cex.main = 1, axes = FALSE)
  box()
  lines(t, yt, col = "darkred", lty=3, lwd=2)
  lines(t, cur_pred, col = "black", lwd=2)
  gen_err[ind] <- mean((cur_pred - yt)^2)
}
```

```{r exportPlot5, echo=FALSE, results='hide', fig.show='hide'}
par(mfrow = c(2,5))
for (ind in seq_along(reg_list)) {
  par(mar = rep(3,4))
  cur_pred <- predict(res[[ind]], as.matrix(t, ncol = 1))
  plot(x, y, pch = "+", lwd = 2, main = print_reg[ind], 
       ylim = range(c(y, yt, cur_pred)), cex.main = 2, axes = FALSE)
  box()
  lines(t, yt, col = "darkred", lty=3, lwd=2)
  lines(t, cur_pred, col = "black", lwd=2)
}
dev.print(png, file = "../figures/simulated_est_decay.png", width = 1200, 
          height = 800, res = 150)
dev.print(pdf, file = "../figures/simulated_est_decay.pdf", width = 12, 
          height = 8)
```

Le graphique ci-dessous illustre précisément l'évolution des erreurs 
d'apprentissage et de test selon la valeur de $-\log_{10}(\lambda)$ (pour 
$\lambda > 0$) et montre qu'un minimum est atteint pour $\lambda=$ 
`r reg_list[which.min(gen_err)]` :

```{r plotErrDecay}
plot(-log10(reg_list[-10]), train_err[-10], pch = "+", col = "darkred",
     type = "b", ylim = c(0, max(c(train_err, gen_err))), 
     xlab = expression(-log[10](lambda)), ylab = "erreur", cex.axis = 1.5,
     cex.lab = 1.5)
lines(-log10(reg_list[-10]), gen_err[-10], col = "black", lty=2)
points(-log10(reg_list[-10]), gen_err[-10], pch = "+", col = "black")
legend(0.7, 0.3, col = c("darkred", "black"), cex = 1.5, bty = "n",
       legend = c("apprentissage", "test"), lwd = 2, lty=c(1,2))
```

```{r exportPlot6, echo=FALSE, results='hide', fig.show='hide'}
plot(-log10(reg_list[-10]), train_err[-10], pch = "+", col = "darkred",
     type = "b", ylim = c(0, max(c(train_err, gen_err))), 
     xlab = expression(-log[10](lambda)), ylab = "erreur", cex.axis = 1.5,
     cex.lab = 1.5)
lines(-log10(reg_list[-10]), gen_err[-10], col = "black", lty=2)
points(-log10(reg_list[-10]), gen_err[-10], pch = "+", col = "black")
legend(0.7, 0.3, col = c("darkred", "black"), cex = 1.5, bty = "n",
       legend = c("apprentissage", "test"), lwd = 2, lty=c(1,2))
dev.print(png, file = "../figures/erreur_decay.png", width = 800, height = 600,
          res = 150)
dev.print(pdf, file = "../figures/erreur_decay.pdf", width=9.6, height=7.2)
```

La valeur de l'erreur de test obtenue pour cet optimum est donnée ci-dessous. On
constatera, en particulier, qu'elle est inférieur à l'erreur de test obtenue 
pour la valeur optimale de $Q$ sans régularisation :
```{r minErrDecay}
min(gen_err)
```


## Calibration du paramètre de régularisation

De la même manière que la calibration du nombre de neurones sur la couche 
cachée a été réalisée, calibrez le paramètre de régularisation $\lambda$ par
validation croisée en utilisant ```tune.nnet```.

Si vous initialisez la graine aléatoire avec ```set.seed(1612)```, vous devriez
obtenir les résultats suivants :

```{r tuneDecay, echo=FALSE}
set.seed(1612)
res_tune <- tune.nnet(x, y, tune.control(nrepeat = iter), size = k,
                      decay = reg_list, linout = TRUE)
```

```{r plotTuneDecay, echo=FALSE}
plot(-log10(reg_list[-10]), res_tune$performances$error[-10], type="b", 
     xlab = expression(-log[10](lambda)), main = "", ylab = "erreur", 
     cex.axis = 1.5, cex.lab = 1.5)
```

```{r exportPlot7, echo=FALSE, results='hide', fig.show='hide'}
plot(-log10(reg_list[-10]), res_tune$performances$error[-10], type="b", 
     xlab = expression(-log[10](lambda)), main = "", ylab = "erreur", 
     cex.axis = 1.5, cex.lab = 1.5)
dev.print(png, file = "../figures/tune_decay.png", width = 800, height = 600, 
          res = 150)
dev.print(pdf, file = "../figures/tune_decay.pdf", width=9.6, height=7.2)
```

```{r bestNNETDecay}
res_tune$best.parameters
res_tune$best.performance
res_tune$best.model
```


# Un exemple sur des données réelles

Enfin, cette dernière partie vous propose de mettre en œuvre ce que vous avez 
appris sur les données ```BostonHousing``` du package mlbench. Ces données 
peuvent être chargées avec les commandes :
```{r loadBoston}
data("BostonHousing")
```
et on trouve des informations sur leur provenance et les questions qu'elles 
soulèvent dans l'aide :
```{r bostonHelp, eval=FALSE}
?BostonHousing
```

En particulier, on cherchera à prédire le revenu médian des habitants des foyers
d'un quartier (```medv```) à partir de caractéristiques socio-économiques de ce 
quartier. Pour simplifier le problème, nous retirons la variable explicative non 
quantitative (```rad```) et effectuons une standardisation des données :
```{r prepareBoston}
df <- BostonHousing
df$chas <- NULL
df <- as.data.frame(scale(df))
```

Utilisez la fonction ```tune.nnet``` pour calibrer simultanément le nombre de
neurones sur la couche cachée parmi $\{5,\, 10,\, 15,\, 20\}$ et le paramètre de
régularisation parmi $\{1,\, 10^{-1},\, 10^{-2},\, 10^{-3}\}$. En fixant la 
graine aléatoire à la valeur ```set.seed(1409)```, vous obtiendrez les résultats
suivants (pour une valeur de ```maxit``` égale à 1000 et 5 répétitions de 
l'apprentissage : **attention !** l'apprentissage est long...) :

```{r bostonNNET, cache=TRUE, echo=FALSE}
set.seed(1409)
res_tune <- tune.nnet(medv ~ ., data = df, size = c(5, 10, 15, 20, 30), 
                      decay = 10^(- (0:3)), linout = TRUE,
                      tunecontrol = tune.control(nrepeat = 5),
                      maxit = 1000)
```

```{r analyseTuneBoston}
plot(res_tune, main = "", color.palette = topo.colors, transform.y = log10)
```

```{r bestBoston}
res_tune$best.parameters
res_tune$best.model
res_tune$best.model$convergence
```

La qualité de l'ajustement peut être visualisée dans le graphique suivant :

```{r predictedBoston}
plot(df$medv, res_tune$best.model$fitted.values, pch = 19, col = "darkred",
     cex = 0.7, xlab = "valeurs réelles", ylab = "valeurs prédites")
abline(0, 1)
```

# Information sur la session **R**

```{r sessionInfo}
session_info()
```
