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.nathalievilla.org/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() ```