library(ggplot2)
library(tidyr)
library(dplyr)
set.seed(123)
n <- 10000
p <- 5
theta_true <- c(-2,-1,0,1,2)

Objectifs du TP

L’objectif de ce TP est d’étudier le comportement d’un estimateur en ligne basé sur le gradient dans les cadres :

On s’intéressera à la convergence de l’estimateur \(\theta_i\) vers le vrai paramètre \(\theta\) et à l’évolution de l’erreur quadratique.

1. Régression linéaire

On considère le modèle

\[ Y_i = X_i^T \theta + \varepsilon_i, \]

avec :

On note \(n = 10000\).

  1. Générer un échantillon \((X_i, Y_i)_{i=1,\dots,n}\).
X <- matrix(rnorm(n*p), n, p)
Y <- X %*% theta_true + rnorm(n)
  1. Écrire une fonction R qui ressorte l’ensemble des estimateurs \(\theta_i\), pour \(i = 0, \ldots, n-1\), définis par la descente de gradient stochastique : \[ \theta_{i+1} = \theta_{i} + \gamma_{i+1} \nabla_\theta \left( (Y_{i+1} - X_{i+1}^T\theta_{i}\right)X_{i+1}. \]
sgd_linear <- function(X, Y, c=1,alpha=0.75) {
  n <- nrow(X); p <- ncol(X)
  theta <- matrix(0, n+1, p)
  for(i in 1:n){
    grad <- -(Y[i] - sum(X[i,]*theta[i,])) * X[i,]
    theta[i+1,] <- theta[i,] - c* i^(-alpha)*grad
  }
  theta
}
theta_path <- sgd_linear(X,Y)
  1. Tracer l’évolution de l’erreur quadratique \(\|\theta_i-\theta\|^2\) pour un échantillon.
err <- rowSums((theta_path - matrix(theta_true,n+1,p,byrow=TRUE))^2)

df <- data.frame(
  i = seq_along(err),
  err = err
)

ggplot(df, aes(x = i, y = err)) +
  geom_line(color = "blue", linewidth = 1) +
  scale_x_log10() +
  scale_y_log10() +
  labs(
    x = "n",
    y = "Erreur quadratique"
  ) +
  theme_minimal()

  1. Répéter l’expérience sur 50 échantillons indépendants et tracer l’erreur quadratique moyenne.
B <- 50
err_mat <- matrix(0,n+1,B)
for(b in 1:B){
  X <- matrix(rnorm(n*p), n, p)
  Y <- X %*% theta_true + rnorm(n)
  th <- sgd_linear(X,Y)
  err_mat[,b] <- rowSums((th - matrix(theta_true,n+1,p,byrow=TRUE))^2)
}
df_mean <- data.frame(
  i = seq_len(nrow(err_mat)),
  err_mean = rowMeans(err_mat)
)

ggplot(df_mean, aes(x = i, y = err_mean)) +
  geom_line(color = "blue", linewidth = 1) +
  scale_x_log10() +
  scale_y_log10() +
  labs(
    x = "n",
    y = "Erreur quadratique moyenne"
  ) +
  theme_minimal()

  1. Refaire la question 4 mais en prenant cette fois-ci \(X \sim \mathcal{N}(0,\text{diag}(1,4,9,16,25))\) et \(c_{\gamma}=0.3,0.1,0.01\).
B <- 50
err_mat_c1 <- matrix(0,n+1,B)
err_mat_c01 <- matrix(0,n+1,B)
err_mat_c001 <- matrix(0,n+1,B)
for(b in 1:B){
  X <- matrix(rnorm(n*p), n, p)%*%diag(1:5)
  Y <- X %*% theta_true + rnorm(n)
  th <- sgd_linear(X,Y,c=0.3)
  th001 <- sgd_linear(X,Y,c=0.01)
  th01 <- sgd_linear(X,Y,c=0.1)
  err_mat_c1[,b] <- rowSums((th - matrix(theta_true,n+1,p,byrow=TRUE))^2)
  err_mat_c01[,b] <- rowSums((th01 - matrix(theta_true,n+1,p,byrow=TRUE))^2)
  err_mat_c001[,b] <- rowSums((th001 - matrix(theta_true,n+1,p,byrow=TRUE))^2)
}
df <- data.frame(
  i = seq_len(nrow(err_mat_c1)),
  `c == 0.3`   =  rowMeans(err_mat_c1),
  `c == 0.1` =  rowMeans(err_mat_c01),
  `c == 0.01`=  rowMeans(err_mat_c001)
)

# Passage en format long
df_long <- df %>%
  pivot_longer(
    cols = -i,
    names_to = "c",
    values_to = "err"
  )
ggplot(df_long, aes(x = i, y = err, color = c)) +
  geom_line(linewidth = 1) +
  scale_y_log10() +
  labs(
    x = "n",
    y = "Erreur quadratique moyenne",
    color = ""
  ) +
  theme_minimal(base_size = 14)

2. Régression logistique

On se place dans le cadre de la régression logistique, i.e. \[ \mathbb{P}(Y_i=1|X_i)=\sigma(X_i^T\theta), \quad \sigma(t)=\frac{1}{1+e^{-t}}. \]

On rappelle que l’algorithme de gradient stochastique est alors défini pour tout \(i=0, \ldots ,n-1\) par \[ \theta_{i+1} = \theta_{i} + \gamma_{t+1} \left( Y_{i+1} - \sigma \left( X_{i+1}^{T}\theta_{i} \right)\right) X_{i+1} \]

Le code suivante permet de générer un échantillon de taille \(n\) pour un \(\theta\) donné.

sigmoid <- function(t) 1/(1+exp(-t))

X <- matrix(rnorm(n*p), n, p)
prob <- sigmoid(X %*% theta_true)
Y <- rbinom(n,1,prob)
  1. Écrire une fonction R qui ressorte l’ensemble des estimateurs
sgd_logistic <- function(X,Y,c=1,alpha=0.75){
  n <- nrow(X); p <- ncol(X)
  theta <- matrix(0,n+1,p)
  for(i in 1:n){
    p_i <- sigmoid(sum(X[i,]*theta[i,]))
    grad <- (p_i - Y[i]) * X[i,]
    theta[i+1,] <- theta[i,] - c*i^(-alpha)*grad
  }
  theta
}
  1. Tracer l’évolution de l’erreur quadratique \(\|\theta_i-\theta\|^2\) pour un échantillon.
theta_path <- sgd_logistic(X,Y)
err <- rowSums((theta_path - matrix(theta_true,n+1,p,byrow=TRUE))^2)
df <- data.frame(
  i = seq_along(err),
  err = err
)

ggplot(df, aes(x = i, y = err)) +
  geom_line(color = "blue", linewidth = 1) +
  scale_y_log10() +
  scale_x_log10() +
  labs(
    x = "n",
    y = "Quadratic error"
  ) +
  theme_minimal()

3. Données MNIST

Introduction aux données MNIST

La base de données MNIST (Modified National Institute of Standards and Technology) est une base de référence très utilisée en apprentissage automatique et en reconnaissance de formes. Elle est composée d’images de chiffres manuscrits allant de 0 à 9 et sert fréquemment de jeu de données “standard” pour tester et comparer des méthodes de classification.

La base MNIST contient : - 60 000 images d’apprentissage - 10 000 images de test

Chaque observation correspond à une image en niveaux de gris de taille 28 × 28 pixels. Les images sont généralement représentées sous forme vectorielle : chaque image est transformée en un vecteur de 784 variables, correspondant aux intensités des pixels, avec des valeurs comprises entre 0 et 255.

L’objectif principal est de construire un modèle capable de prédire le chiffre manuscrit représenté sur une image à partir des valeurs des pixels. Il s’agit donc d’un problème de classification supervisée, le plus souvent : - multiclasse, avec 10 classes (chiffres de 0 à 9), - ou parfois binaire, par exemple pour distinguer un chiffre donné des autres.

Ici, on se placera dans le cadre binaire, i.e. on ne considèrera ques les chiffres 0 et 1 et on appliquera la régression logistique pour faire de la prévision.

Le code suivant permet d’extraire les données d’entrainement :

library(dslabs)
mnist <- read_mnist()


data_train=mnist$train
I=which(data_train$label == 0 | data_train$label == 1)

Y_train=data_train$labels[I]
X_train=data_train$images[I,]

On obtient ce type d’images pour les données MNIST :

par(mfrow = c(1, 2), mar = c(1, 1, 2, 1))

img1 <- matrix(X_train[1, ], 28, 28, byrow = TRUE)
img2 <- matrix(X_train[2, ], 28, 28, byrow = TRUE)

image(
  img1[, 28:1],
  col = gray.colors(256),
  axes = FALSE,
  asp = 1,
  main = paste("Label =", Y_train[1])
)

image(
  img2[, 28:1],
  col = gray.colors(256),
  axes = FALSE,
  asp = 1,
  main = paste("Label =", Y_train[2])
)

Dans ce qui suit, on se place dans le cadre de la régression logistique, i.e. on suppose \[ Y_{i} | X_{i} \sim \mathcal{B}\left(\sigma \left( \theta^{T}X_{i} \right) \right) . \]

  1. Estimer le paramètre \(\theta\) à l’aide de l’algorithme de gradient stochastique.
train_result=sgd_logistic(X_train,Y_train)
  1. Pour prédire \(Y_{i}\), on considère l’estimateur \(\hat{\theta}\), et la prévision \(\hat{Y}_{i} = 1\) si \(\sigma \left( X_{i}^{T}\hat{\theta}\right) \geq 1/2\) et \(0\) sinon. Calculer vos prévision sur les données d’entraînement, puis visualiser l’erreur de prévision avec la fonction table avant de calculer le taux de bonnes prévisions.
theta_final=train_result[length(Y_train),]
Y_prev=rep(0,length(Y_train))
for (i in 1:length(Y_train))
{
  sig=sigmoid(sum(theta_final*X_train[i,]))
  if (sig > 0.5)
  {
    Y_prev[i]=1
  }
}

tab=table(Y_train,Y_prev)
score=sum(diag(tab))/length(Y_train)

tab
##        Y_prev
## Y_train    0    1
##       0 5887   36
##       1   42 6700
score
## [1] 0.9938413
  1. Le code suivant permet d’extraire les données test.
data_test=mnist$test
J=which(data_test$label == 0 | data_test$label == 1)

Y_test=data_test$labels[J]
X_test=data_test$images[J,]

Calculer les prédictions pour les données tests. Attention, il ne faut pas rééstimer \(\theta\) mais conserver l’estimateur précédent.

Y_prev_test=rep(0,length(Y_test))
for (i in 1:length(Y_test))
{
  sig=sigmoid(sum(theta_final*X_test[i,]))
  if (sig > 0.5)
  {
    Y_prev_test[i]=1
  }
}


tab_test=table(Y_test,Y_prev_test)

score_test=sum(diag(tab_test))/length(Y_test)

tab_test
##       Y_prev_test
## Y_test    0    1
##      0  976    4
##      1    3 1132
score_test
## [1] 0.9966903