Parcours-ML/05_Metamethods/exemples.Rmd

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