|
|
@@ -0,0 +1,189 @@ |
|
|
|
```{r libs} |
|
|
|
library(tidyverse) |
|
|
|
``` |
|
|
|
|
|
|
|
# Cross-validation |
|
|
|
|
|
|
|
## Création du set de données |
|
|
|
|
|
|
|
On initialise le RNG avec une graine fixe. |
|
|
|
On génère une distribution uniforme sur l'axe des **x**, |
|
|
|
et **y** selon un polynôme de degré 2 avec un bruit gaussien. |
|
|
|
|
|
|
|
```{r dataset CV} |
|
|
|
set.seed(1234) |
|
|
|
tibble(x = runif(100, 0, 20), |
|
|
|
y = x ^ 2 + 2 * x + rnorm(100, 0, 4)) -> df |
|
|
|
``` |
|
|
|
|
|
|
|
## Partition train/test |
|
|
|
|
|
|
|
Partition du set de données en un set d'**entraînement** (80%), |
|
|
|
et un set de **test** (les 20% restants). |
|
|
|
|
|
|
|
```{r partition CV} |
|
|
|
df %>% |
|
|
|
sample_n(80) -> df_train |
|
|
|
|
|
|
|
df %>% |
|
|
|
anti_join(df_train) -> df_test |
|
|
|
``` |
|
|
|
|
|
|
|
## Fonctions |
|
|
|
|
|
|
|
Création de quelques fonctions pour la réalisation de la cross-validation. |
|
|
|
|
|
|
|
### Générateur de modèles |
|
|
|
|
|
|
|
Fonction prenant une dataframe et un degré de polynôme pour générer le modèle correspondant. |
|
|
|
La syntaxe `I(...)` dans un modèle sert à décrire des calculs littéraux sur les variables. |
|
|
|
Ici, on concatène `x^n` avec `n` allant de 1 au degré du polynôme, et on applique le modèle. |
|
|
|
|
|
|
|
```{r model function} |
|
|
|
model_lm <- function(degre) |
|
|
|
{ |
|
|
|
polynome <- str_c("I(x^", 1:degre, ")", collapse = " + ") |
|
|
|
formule <- str_c("y ~ ", polynome) |
|
|
|
|
|
|
|
# lm(as.formula(formule), data = df) |
|
|
|
partial(lm, formula = as.formula(formule)) |
|
|
|
} |
|
|
|
``` |
|
|
|
|
|
|
|
### Validateur de modèles |
|
|
|
|
|
|
|
Fonction prenant un set d'entraînement, un set de test, et un degré de polynôme, et calculant la _Mean Square Error_ du modèle polynomial correspondant. |
|
|
|
|
|
|
|
```{r MSE of a model} |
|
|
|
mse <- function(train, test, model) |
|
|
|
{ |
|
|
|
modele <- model(train) |
|
|
|
|
|
|
|
verite <- test$y |
|
|
|
predictions <- predict(modele, test) |
|
|
|
|
|
|
|
mean((verite - predictions) ^ 2) |
|
|
|
} |
|
|
|
``` |
|
|
|
|
|
|
|
### Subsampling pour CV |
|
|
|
|
|
|
|
Fonction prenant un dataset et un nombre **k** comme arguments, et générant **k** permutations train/test. |
|
|
|
On passe par la création d'une variable _groupVar_ qui est un tirage aléatoire équilibré (le vecteur 1:k est répété sur la longueur du dataset, puis mélangé). |
|
|
|
Ensuite, chaque k-ième du dataset est pris comme set de test, pendant que le reste sert de set d'entraînement. |
|
|
|
|
|
|
|
```{r subsampling} |
|
|
|
subsampl <- function(df, k) |
|
|
|
{ |
|
|
|
df %>% |
|
|
|
mutate(groupVar = rep(1:k, |
|
|
|
length.out = nrow(df)) %>% |
|
|
|
sample) %>% |
|
|
|
group_nest(groupVar, .key = "test") %>% |
|
|
|
mutate(train = test %>% map(~ anti_join(df, .))) %>% |
|
|
|
select(-groupVar) |
|
|
|
} |
|
|
|
``` |
|
|
|
|
|
|
|
### Cross-validation d'un modèle |
|
|
|
|
|
|
|
Fonction appliquant la procédure de cross-validation pour un modèle. |
|
|
|
On créé tout d'abord les k permutations, puis on calcule la **MSE** pour chaque permutation, avec le modèle donné. |
|
|
|
Enfin, on calcule la moyenne des MSE de chacune des permutations. |
|
|
|
|
|
|
|
```{r CV_one} |
|
|
|
CV_one <- function(df, model, k) |
|
|
|
{ |
|
|
|
df %>% |
|
|
|
subsampl(k) %>% |
|
|
|
mutate(mse = map2_dbl(train, test, ~mse(.x, .y, model))) %>% |
|
|
|
pull(mse) %>% |
|
|
|
mean |
|
|
|
} |
|
|
|
``` |
|
|
|
|
|
|
|
### Cross-validation sur tous les modèles |
|
|
|
|
|
|
|
Fonction appliquant la procédure de cross-validation à tous les modèles en argument. |
|
|
|
Le dataset est dupliqué pour chaque modèle, et on applique ensuite la procédure unique pour chaque couple degré-modèle. |
|
|
|
|
|
|
|
```{r CV} |
|
|
|
CV <- function(df, k, models) |
|
|
|
{ |
|
|
|
df %>% |
|
|
|
crossing(degre = 1:length(models)) %>% |
|
|
|
group_nest(degre) %>% |
|
|
|
mutate(mse = map2_dbl(data, degre, ~ CV_one(.x, .y, k))) |
|
|
|
} |
|
|
|
``` |
|
|
|
|
|
|
|
## Application |
|
|
|
|
|
|
|
Il ne reste plus qu'à appliquer la cross-validation souhaitée sur notre set d'entraînement, et choisir le modèle le plus approprié. |
|
|
|
|
|
|
|
```{r application} |
|
|
|
2:5 %>% |
|
|
|
map(model_lm) -> liste_modeles |
|
|
|
|
|
|
|
df_train %>% |
|
|
|
CV(80, liste_modeles) |
|
|
|
``` |
|
|
|
|
|
|
|
## Modèle final |
|
|
|
|
|
|
|
```{r modèle final} |
|
|
|
lm(y ~ x, data = df_train) -> fit |
|
|
|
|
|
|
|
predict(fit, df_test) -> predictions |
|
|
|
|
|
|
|
mean((df_test$y - predictions)^2) |
|
|
|
``` |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
# Bootstrap |
|
|
|
|
|
|
|
Le bootstrap est une méthode de _resampling_ qui consiste à tirer un grand nombre **n** de tirages au sort avec remise de la distribution en entrée. |
|
|
|
|
|
|
|
## Bootstrap de la moyenne |
|
|
|
|
|
|
|
Fonction prenant un vecteur à bootstrap, et le nombre de répétitions. |
|
|
|
Il s'agit d'un simple tirage avec remise, effectué **b** fois, dont on calcule la moyenne. |
|
|
|
|
|
|
|
```{r bootstrap} |
|
|
|
bootstrap <- function(x, n) |
|
|
|
{ |
|
|
|
replicate(n, sample(x, replace = T) %>% mean) |
|
|
|
} |
|
|
|
``` |
|
|
|
|
|
|
|
## Exemple |
|
|
|
|
|
|
|
Calcul de l'intervalle à 95% de la moyenne. |
|
|
|
|
|
|
|
### Données |
|
|
|
|
|
|
|
Distribution normale aléatoire de 30 sujets. |
|
|
|
|
|
|
|
```{r} |
|
|
|
x <- rnorm(30) |
|
|
|
qplot(x) |
|
|
|
mean(x) |
|
|
|
t.test(x)$conf.int |
|
|
|
``` |
|
|
|
|
|
|
|
### Production du bootstrap |
|
|
|
|
|
|
|
```{r apply bootstrap} |
|
|
|
b <- bootstrap(x, 10000) |
|
|
|
|
|
|
|
tibble(x = b) %>% |
|
|
|
ggplot() + |
|
|
|
aes(x = x) + |
|
|
|
geom_histogram() + |
|
|
|
geom_vline(xintercept = mean(b)) + |
|
|
|
geom_vline(xintercept = quantile(b, probs = c(.025, .975))) |
|
|
|
|
|
|
|
mean(b) |
|
|
|
quantile(b, probs = c(.025, .975)) |
|
|
|
``` |