diff --git a/21_b_tp_projection_aleatoire_housing_pca.Rmd b/21_b_tp_projection_aleatoire_housing_pca.Rmd new file mode 100644 index 0000000..53f9bc9 --- /dev/null +++ b/21_b_tp_projection_aleatoire_housing_pca.Rmd @@ -0,0 +1,237 @@ +# TP - Projection aléatoire - `housing` PCA - Correction + +```{r, include=FALSE} +source("05_svd_pca.R", local = knitr::knit_global()) +``` + +```{r} +set.seed(1123) +``` + +## Présentation du jeu de données `housing` + +`housing` est un jeu de données célèbre aux nombreuses vertues pédagogiques^[https://www.kaggle.com/datasets/harrywang/housing]. Il permet d'expérimenter sur un problème de régression réaliste, viz. prédire la valeur médiane d'une maison en fonction des caractéristiques de son quartier. Après une phase d'exploration des données, comparer, sur un jeu de test, un modèle linéaire par régression ridge et un modèle par projection aléatoire. + +```{r cache=TRUE} +housing <- read.csv(file="data/housing.csv",header=TRUE) +str(housing) +``` + +Variable | Type | Commentaire +-------------------|---------|------------------------------------ +longitude | numeric | élevée pour un quartier à l'ouest +latitude | numeric | élevée pour un quartier au nord +housing_median_age | numeric | âge médian d'une maison du quartier +total_rooms | numeric | total des pièces pour les maisons du quartier +total_bedrooms | numeric | total des chambres pour les maisons du quartier +population | numeric | nombre de résidents du quartier +households | numeric | nombre de familles du quartier +median_income | numeric | revenu médian des familles du quartier (USD 10k) +median_house_value | numeric | valeur médiane d'une maison du quartier (USD) +ocean_proximity | factor | estimation de la distance à l'océan + +## Prétraitements + +Nous transformons la variable `ocean_proximity`, pour l'instant représentée par une chaîne de caractères, en une variable catégorielle (ou facteur) qui sera représentée par un codage disjonctif complet (ou _one hot encoding_). + +```{r} +housing$ocean_proximity <- as.factor(housing$ocean_proximity) +summary(housing) +``` + +Nous remplaçons les valeurs manquantes pour la variable `total_bedrooms` par la médiane. + +```{r} +housing$total_bedrooms[is.na(housing$total_bedrooms)] <- median(housing$total_bedrooms, na.rm=TRUE) +sum(is.na(housing)) +``` + +La catégorie `ISLAND` de la variable `ocean_proximity` est extrêmement sous représentée. Nous proposons de supprimer les observations qui lui correspondent. + +```{r} +housing <- housing[housing$ocean_proximity != "ISLAND", ] +nrow(housing) +``` + +## Analyse exploratoire par PCA + +Explorons les variables numériques par analyse en composantes principales. + +Nous sélectionnons les variables explicatives numériques dans une matrice $\mathbf{X}$. +```{r} +X <- housing[,c("longitude", "latitude", "housing_median_age", "total_rooms", + "total_bedrooms", "population", "households", "median_income", + "median_house_value")] +``` + +### Première analyse + +Nous réalisons une première fois l'analyse en composantes principales. +```{r} +fam <- fa(X) # fam pour 'factor analysis model' +``` + +Nous affichons le pourcentage de variance expliquée par chaque axe principal. +```{r} +fam$prctPrcp +``` + +Nous affichons, sur les deux premiers axes principaux, les centres des clusters qui contribuent le plus à ces axes. +```{r} +print(fam) +``` + +```{r} +far <- away(fam) +``` + +Le cluster `r far$id` (`far$id`) contribue à expliquer `r round(fam$ctr[far$id,1]*100, 2)` % (`round(fam$ctr[far$id,1]*100, 2)`) de la variance du premier axe. Il est composé de `r far$size` (`far$size`) quartiers : +```{r} +housing[far$names,] +``` + +Ces quartiers semblent remarquables par le grand nombre de familles qui les habitent. +```{r} +boxplot(housing$households) +``` + +Nous remarquons que le cluster `r far$id` est principalement expliqué par le premier axe principal. +```{r} +fam$cos2[far$id,] +``` + +Les contributions des variables aux axes principaux nous montrent que le premier axe principal correspond surtout à un facteur taille de variables très corrélées : `total_rooms`, `total_bedrooms`, `population` et `households`. +```{r} +fam$varctr +``` + +### Seconde analyse + +Pour mieux visualiser des phénomènes intéressants qui seraient sinon masqués par ce facteur taille, nous proposons d'introduire de nouvelles variables pour décrire le nombre de pièces, le nombre de chambres et le nombre de personnes relativement au nombre de familles du quartier. +```{r} +housing['rooms_per_household'] <- housing['total_rooms'] / housing['households'] +housing['bedrooms_per_household'] <- housing['total_bedrooms'] / housing['households'] +housing['population_per_household'] <- housing['population'] / housing['households'] +``` + +Nous recommençons l'analyse avec ces nouvelles variables. +```{r} +X <- housing[,c("longitude", "latitude", "housing_median_age", "rooms_per_household", + "bedrooms_per_household", "population_per_household", "median_income", + "median_house_value")] +fam <- fa(X) +``` + +Nous affichons le pourcentage de variance expliquée par chaque axe principal. +```{r} +fam$prctPrcp +``` + +Nous affichons, sur les deux premiers axes principaux, les centres des clusters qui contribuent le plus à ces axes. +```{r} +print(fam) +``` + +```{r} +far <- away(fam) +``` + +Le cluster `r far$id` contribue à expliquer `r round(fam$ctr[far$id,1]*100, 2)` % de la variance du premier axe. Il est composé de `r far$size` quartiers : +```{r} +housing[far$names,] +``` + +Cependant, la contribution du cluster `r far$id` à la variance est surtout expliquée par l'axe $3$. Ce dernier est lui-même surtout expliqué par les variables `rooms_per_household` et `bedrooms_per_household` (elles-mêmes très corrélées, ce n'est pas étonnant). +```{r} +fam$cos2[far$id,] +``` + +```{r} +fam$varctr +``` + +Nous vérifions que ces quartiers correspondent effectivement à des valeurs très atypiques de la variable `rooms_per_household`. +```{r} +boxplot(housing$rooms_per_household) +``` + +Nous remarquons aussi que la valeurs des maisons sur l'un de ces quartiers est très élevée. En fait, c'est la plus grande valeur rencontrée sur ce jeu de données. +```{r} +capMedianHouseValue <- 500001 +nbCapMedianHouseValue <- table(housing$median_house_value)[as.character(capMedianHouseValue)] +``` + +Pour `r nbCapMedianHouseValue` quartiers, la valeur de la variable `median_house_value` est égale à `r capMedianHouseValue`. Il semble que les maisons dont la valeur dépasse une certaine somme aient été toutes enregistrées à cette somme maximale. Ces quartiers risquent d'avoir une mauvaise influence sur notre modèle prédictif. Nous proposons donc de les retirer. + +### Troisième analyse + +```{r} +housing <- housing[housing$median_house_value < capMedianHouseValue, ] +nrow(housing) +``` + +```{r} +X <- housing[,c("longitude", "latitude", "housing_median_age", "rooms_per_household", + "bedrooms_per_household", "population_per_household", "median_income", + "median_house_value")] +fam <- fa(X) +``` + +Nous affichons le pourcentage de variance expliquée par chaque axe principal. +```{r} +fam$prctPrcp +``` + +Nous vérifions les contributions des variables aux axes principaux. +```{r} +fam$varctr +``` + +Nous affichons, sur les deux premiers axes principaux, les centres des clusters qui contribuent le plus à ces axes. +```{r} +print(fam) +``` + +```{r} +far <- away(fam) +``` + +Le cluster `r far$id` contribue à expliquer `r round(fam$ctr[far$id,1]*100, 2)` % de la variance du premier axe. Il est composé de `r far$size` quartiers : +```{r} +housing[far$names,] +``` + +Ce quartier, que nous avons déjà rencontré dans la précédente analyse est vraiment atypique. Il s'agit d'un petit quartier avec énormément de pièces par foyer. Nous proposons de le retirer. +```{r} +boxplot(housing$bedrooms_per_household) +``` + +### Quatrième analyse + +```{r} +housing <- housing[!(row.names(housing) %in% far$names),] +``` + +```{r} +X <- housing[,c("longitude", "latitude", "housing_median_age", "rooms_per_household", + "bedrooms_per_household", "population_per_household", "median_income", + "median_house_value")] +fam <- fa(X) +``` + +Nous affichons le pourcentage de variance expliquée par chaque axe principal. +```{r} +fam$prctPrcp +``` + +Nous vérifions les contributions des variables aux axes principaux. +```{r} +fam$varctr +``` + +Nous affichons, sur les deux premiers axes principaux, les centres des clusters qui contribuent le plus à ces axes. +```{r} +print(fam) +``` + +Nous proposons d'arrêter ici l'analyse exploratoire. \ No newline at end of file diff --git a/21_c_tp_projection_aleatoire_housing_elm.Rmd b/21_c_tp_projection_aleatoire_housing_elm.Rmd new file mode 100644 index 0000000..85a6800 --- /dev/null +++ b/21_c_tp_projection_aleatoire_housing_elm.Rmd @@ -0,0 +1,144 @@ +# TP - Projection aléatoire - `housing` ELM - Correction + +```{r, include=FALSE} +source("15_loocv.R", local = knitr::knit_global()) +``` + +```{r} +set.seed(1123) +``` + +## Créations des jeux d'entraînement et de test + +Nous reproduisons les prétraitements découverts pendant l'analyse exploratoire +```{r} +housing <- read.csv(file="data/housing.csv",header=TRUE) +housing$ocean_proximity <- as.factor(housing$ocean_proximity) +housing$total_bedrooms[is.na(housing$total_bedrooms)] <- + median(housing$total_bedrooms, na.rm=TRUE) +housing <- housing[housing$ocean_proximity != "ISLAND", ] +housing['rooms_per_household'] <- housing['total_rooms'] / housing['households'] +housing['bedrooms_per_household'] <- housing['total_bedrooms'] / housing['households'] +housing['population_per_household'] <- housing['population'] / housing['households'] +housing <- housing[housing$median_house_value < 500001, ] +housing <- housing[-1980,] +``` + +Nous rappelons une fonction pour créer un jeu d'entraînement et un jeu de test. +```{r} +# Séparer le jeu de données en un jeu d'entraînement et un jeu de test +# INPUT : jeu de données initial et proportion des données conservées pour +# l'entraînement. +splitdata <- +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]), + test = list(X = data$X[-entridx,,drop=FALSE], Y = data$Y[-entridx])) +} +``` + +Nous créons les jeux de données. +```{r} +data <- list(X=housing[,c("longitude", "latitude", "housing_median_age", + "rooms_per_household", "bedrooms_per_household", + "population_per_household", "households", + "median_income")], + Y=housing[,c("median_house_value")]) +splitres <- splitdata(data,0.8) +entr <- splitres$entr +test <- splitres$test +``` + +Si nous essayons d'inférer un modèle ELM avec l'appel ci-dessous qui utilise le code développé pour les données synthétiques `sinc`, nous rencontrons un bug car la projection aléatoire des données peut créer des colonnes avec seulement des valeurs négative. Lorsqu'une telle colonne est transformée, par exemple par la fonction non linéaire `ReLU`, elle deviendra pleine de zéros et rendra impossible le calcul de la décomposition en valeurs singulières. C'est pourquoi, pour les modèles qui emploient des fonctions de transformation non linéaires, les données sont souvent normalisées entre $0$ et $1$. Par ailleurs, s'il reste des colonnes nulles dans $\mathbf{H}$, nous les retirons (et nous mettons à jour en conséquence les vecteurs aléatoires $\mathbf{W}$ et $\mathbf{b}$). +```{r eval=FALSE} +rm <- elm(entr$X,entr$Y) +``` + +## Modèle ELM + +```{r} +relu <- function(x){ifelse(x>=0,x,0)} +logistic <- function(x){1/(1+exp(-x))} + +elm <- +function(X,y,m=1000,method=c('logistic', 'ReLU'),lambdas=NULL) +{ + method <- match.arg(method) + + X <- as.matrix(X) + Xisnum <- apply(X,2,is.numeric) # les colonnes qui sont des facteurs + # n'ont pas besoin d'être normalisées + Xmax <- apply(X[,Xisnum],2,max) + Xmin <- apply(X[,Xisnum],2,min) + Xrange <- Xmax - Xmin + X[,Xisnum] <- sweep(X[,Xisnum],2,Xmin,'-') + X[,Xisnum] <- sweep(X[,Xisnum],2,Xrange,'/') + + p <- ncol(X) + W <- matrix(runif(p*m, min=-1, max=1), nrow=p, ncol=m) + H <- X%*%W + b <- runif(m) + H <- sweep(H,2,b,'+') + + if(method=='logistic') { + H <- logistic(H) + } else if(method=='ReLU') { + H <- relu(H) + } + + nullcols <- which(apply(H,2,sum)<=1) + if(length(nullcols)>0) { + H <- H[,-nullcols] + W <- W[,-nullcols] + b <- b[-nullcols] + m <- m - length(nullcols) + } + + rm <- ridge(H,y,lambdas) + r <- list(ridge = rm, + m = m, + W = W, + b = b, + method=method, + isnum = Xisnum, + min = Xmin, + range = Xrange) + class(r) <- "elm" + return(r) +} + +predict.elm <- +function(o, newdata) +{ + newdata <- as.matrix(newdata) + newdata[,o$isnum] <- sweep(newdata[,o$isnum],2,o$min,'-') + newdata[,o$isnum] <- sweep(newdata[,o$isnum],2,o$range,'/') + + H <- newdata %*% o$W + H <- sweep(H,2,o$b,'+') + if(o$method=='logistic') { + H <- logistic(H) + } else if(o$method=='ReLU') { + H <- relu(H) + } + yh <- predict(o$ridge,H) +} +``` + +## Prédictions sur le jeu de test + +Nous apprenons un modèle par projection aléatoire et un modèle ridge. +```{r} +lambdas <- 10^seq(-3,2,length.out=10) +rm <- elm(entr$X,entr$Y,m=200,method='logistic',lambdas=lambdas) +lm <- ridge(entr$X,entr$Y) +rmYh <- predict(rm,test$X) +lmYh <- predict(lm,test$X) +rmTestMAE <- mean(abs(rmYh - test$Y)) +lmTestMAE <- mean(abs(lmYh - test$Y)) +``` + +Sur le jeu de test, le modèle par projection aléatoire commet une erreur absolue moyenne de `r rmTestMAE` tandis que le modèle linéaire régularisé commet une erreur de `r lmTestMAE`. \ No newline at end of file diff --git a/21_tp_projection_aleatoire.Rmd b/21_tp_projection_aleatoire.Rmd new file mode 100644 index 0000000..9093e3c --- /dev/null +++ b/21_tp_projection_aleatoire.Rmd @@ -0,0 +1,126 @@ +# Projection aléatoire - `sinc` - Correction + +```{r, include=FALSE} +source("15_loocv.R", local = knitr::knit_global()) +``` + +```{r} +set.seed(1123) +``` + +L'idée de cette expérimentation provient de la la référence [@huang2006extreme] qui introduit le modèle prédictif dit _Extreme Learning Machine_. + +## Générer un jeu de données synthétique avec la fonction `sinc` + +Utiliser la fonction `sinc` pour générer un jeu de données $\{y_i,x_i\}$. +$$ +y(x)= +\begin{cases} +sin(x)/x & x\neq0, \\ +1 & x=0 +\end{cases} +$$ + +Les $x_i$ sont distribués de façon uniforme sur l'intervalle $(-10,10)$. Un bruit (par exemple uniforme ou gaussien) est ajouté aux étiquettes $y_i$. + +Nous commençons par définir la fonction `sinc`. +```{r} +sinc <- +function(x) +{ + y <- sin(x) / x + y[x==0] <- 1 + return(y) +} +``` + +Puis nous créons le jeu de données. +```{r} +n <- 100 +X <- runif(n, min=-10, max=10) +X <- as.matrix(X) +y <- sinc(X) + runif(n, min=-0.2, max=0.2) +``` + +Et nous l'affichons. +```{r} +xplt <- seq(-10,10,length.out=100) +plot(xplt,sinc(xplt), type='l') +points(X,y, pch=20) +``` + +## Régression ridge après projections non-linéaires aléatoires + +Créer un modèle qui a la forme d'un réseau de neurones à une couche. +$$ +\hat{y}_i = \sum_{j=1}^{m} \beta_j g\left( \mathbf{w_j}\mathbf{x_i}^T + b_j \right) \quad i=1,\dots,n +$$ +Cependant, contrairement à un réseau de neurones, les $m$ vecteurs $\mathbf{w_j}$ et scalaires $b_j$ sont initialisés aléatoirement et ne sont jamais modifiés. La fonction $g$ est une transformation non-linéaire (par ex., la fonction ReLU $g(x)=max(0,x)$, ou la fonction sigmoïde $g(x)=1/(1+exp(-x))$, etc.). Seuls les paramètres $\beta_j$ sont inférés par régression ridge en utilisant la matrice $\mathbf{H}$ ci-dessous au lieu de la martice $\mathbf{X}$ initiale. +$$ +\mathbf{H} = +\left[ \begin{array}{ccc} +g\left( \mathbf{w_1}\mathbf{x_1}^T + b_1 \right) & \dots & g\left( \mathbf{w_m}\mathbf{x_1}^T + b_m \right) \\ +\vdots & \dots & \vdots \\ +g\left( \mathbf{w_1}\mathbf{x_n}^T + b_1 \right) & \dots & g\left( \mathbf{w_m}\mathbf{x_n}^T + b_m \right) \\ +\end{array} \right] +$$ + +Étudier les effets de la taille du jeu d'apprentissage, de la quantité de bruit, de l'hyperparamètre de régularisation, du nombre de transformations aléatoires, de la fonction non linéaire choisie, etc. + +Nous commençons par créer le modèle en étendant la régression ridge. +```{r} +relu <- function(x){ifelse(x>=0,x,0)} +logistic <- function(x){1/(1+exp(-x))} + +elm <- +function(X,y,m=2000,method=c('logistic', 'ReLU'),lambdas=NULL) +{ + method <- match.arg(method) + + X <- as.matrix(X) + p <- ncol(X) + W <- matrix(runif(p*m, min=-1, max=1), nrow=p, ncol=m) + H <- X%*%W + b <- runif(m) + H <- sweep(H,2,b,'+') + + if(method=='logistic') { + H <- logistic(H) + } else if(method=='ReLU') { + H <- relu(H) + } + + rm <- ridge(H,y,lambdas) + r <- list(ridge = rm, + m = m, + method=method, + W = W, + b = b) + class(r) <- "elm" + return(r) +} + +predict.elm <- +function(o, newdata) +{ + newdata <- as.matrix(newdata) + H <- newdata %*% o$W + H <- sweep(H,2,o$b,'+') + if(o$method=='logistic') { + H <- logistic(H) + } else if(o$method=='ReLU') { + H <- relu(H) + } + + yh <- predict(o$ridge,H) +} +``` + +Nous appliquons le modèle au jeu de données synthétique. +```{r} +rm <- elm(X,y,m=500,method='logistic') +yh <- predict(rm,xplt) +plot(xplt,sinc(xplt),type='l',col='blue') +lines(xplt,yh,type='l',col='green') +points(X,y, pch=20) +``` diff --git a/21_tp_projection_aleatoire_sujet.Rmd b/21_tp_projection_aleatoire_sujet.Rmd new file mode 100644 index 0000000..885cc65 --- /dev/null +++ b/21_tp_projection_aleatoire_sujet.Rmd @@ -0,0 +1,60 @@ +# TP - Projection aléatoire - Sujet + +```{r} +set.seed(1123) +``` + +L'idée de cette expérimentation provient de la la référence [@huang2006extreme] qui introduit le modèle prédictif dit _Extreme Learning Machine_. + +## Générer un jeu de données synthétique avec la fonction `sinc` + +Utiliser la fonction `sinc` pour générer un jeu de données $\{y_i,x_i\}$. +$$ +y(x)= +\begin{cases} +sin(x)/x & x\neq0, \\ +1 & x=0 +\end{cases} +$$ + +Les $x_i$ sont distribués de façon uniforme sur l'intervalle $(-10,10)$. Un bruit (par exemple uniforme ou gaussien) est ajouté aux étiquettes $y_i$. + +## Régression ridge après projections non-linéaires aléatoires + +Créer un modèle qui a la forme d'un réseau de neurones à une couche. +$$ +\hat{y}_i = \sum_{j=1}^{m} \beta_j g\left( \mathbf{w_j}\mathbf{x_i}^T + b_j \right) \quad i=1,\dots,n +$$ +Cependant, contrairement à un réseau de neurones, les $m$ vecteurs $\mathbf{w_j}$ et scalaires $b_j$ sont initialisés aléatoirement et ne sont jamais modifiés. La fonction $g$ est une transformation non-linéaire (par ex., la fonction ReLU $g(x)=max(0,x)$, ou la fonction sigmoïde $g(x)=1/(1+exp(-x))$, etc.). Seuls les paramètres $\beta_j$ sont inférés par régression ridge en utilisant la matrice $\mathbf{H}$ ci-dessous au lieu de la martice $\mathbf{X}$ initiale. +$$ +\mathbf{H} = +\left[ \begin{array}{ccc} +g\left( \mathbf{w_1}\mathbf{x_1}^T + b_1 \right) & \dots & g\left( \mathbf{w_m}\mathbf{x_1}^T + b_m \right) \\ +\vdots & \dots & \vdots \\ +g\left( \mathbf{w_1}\mathbf{x_n}^T + b_1 \right) & \dots & g\left( \mathbf{w_m}\mathbf{x_n}^T + b_m \right) \\ +\end{array} \right] +$$ + +Étudier les effets de la taille du jeu d'apprentissage, de la quantité de bruit, de l'hyperparamètre de régularisation, du nombre de transformations aléatoires, de la fonction non linéaire choisie, etc. + +## Jeu de données `housing` + +`housing` est un jeu de données célèbre aux nombreuses vertues pédagogiques^[https://www.kaggle.com/datasets/harrywang/housing]. Il permet d'expérimenter sur un problème de régression réaliste, viz. prédire la valeur médiane d'une maison en fonction des caractéristiques de son quartier. Après une phase d'exploration des données, comparer, sur un jeu de test, un modèle linéaire par régression ridge et un modèle par projection aléatoire. + +```{r cache=TRUE} +housing <- read.csv(file="data/housing.csv",header=TRUE) +str(housing) +``` + +Variable | Type | Commentaire +-------------------|---------|------------------------------------ +longitude | numeric | élevée pour un quartier à l'ouest +latitude | numeric | élevée pour un quartier au nord +housing_median_age | numeric | âge médian d'une maison du quartier +total_rooms | numeric | total des pièces pour les maisons du quartier +total_bedrooms | numeric | total des chambres pour les maisons du quartier +population | numeric | nombre de résidents du quartier +households | numeric | nombre de familles du quartier +median_income | numeric | revenu médian des familles du quartier (USD 10k) +median_house_value | numeric | valeur médiane d'une maison du quartier (USD) +ocean_proximity | factor | évaluation grossière de la distance à l'océan \ No newline at end of file