p6e7p7/intro_to_ml
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

ajout d'un TP pour expérimenter sur la projection aléatoire et l'analyse factorielle apr PCA sur le jeu de données California Housing

Browse Source
This commit is contained in:
Pierre-Edouard Portier 2022-04-16 12:01:04 +02:00
parent b96da84a8d
commit 01c918b1d1
4 changed files with 567 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

237
21_b_tp_projection_aleatoire_housing_pca.Rmd Normal file
View File

@ -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.

144
21_c_tp_projection_aleatoire_housing_elm.Rmd Normal file
View File

@ -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`.

126
21_tp_projection_aleatoire.Rmd Normal file
View File

@ -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)
```

60
21_tp_projection_aleatoire_sujet.Rmd Normal file
View File

@ -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