Cours

Cours.Semaine10TME10 History

Hide minor edits - Show changes to output

Changed line 26 from:
* Vous obtenez un écart-type {$\sigma$} pas simple multiplication sur le tirage précédent
to:
* Vous obtenez un écart-type {$\sigma$} par simple multiplication sur le tirage précédent
Deleted line 19:
c = 1
Changed line 79 from:
## {$w_{t+1} \leftarrow w_{t} - \epsilon \nabla_w C(w)$}
to:
## {$w_{t+1} \leftarrow w_{t} - \epsilon \nabla_w C(w_t)$}
Changed lines 139-140 from:
!!!%subsection% Extension non-linéaire
to:
!!!%subsection% Extension non-linéaire (solution analytique)
Changed line 147 from:
* estimer un vecteur w sur ces données
to:
* estimer un vecteur w sur ces données en utilisant la solution analytique (les données sont petites, pas la peine de s'embéter)
Added lines 130-132:
costPath = np.array([np.log(((X.dot(wtmp)-y)**2).sum()) for wtmp in allw])
costOpt  = np.log(((X.dot(wstar)-y)**2).sum())

Added lines 136-137:
ax.scatter(wstar[0], wstar[1],costOpt, c='r')
ax.plot(allw[:,0],allw[:,1],costPath, 'b+-' ,lw=2 )
Changed lines 84-85 from:
wstar = np.linalg.solve(X.T.dot(X), X.T.dot(ylin)) # pour se rappeler du w optimal
to:
wstar = np.linalg.solve(X.T.dot(X), X.T.dot(y)) # pour se rappeler du w optimal
Changed line 113 from:
cost = np.array([[np.log(((X.dot(np.array([w1i,w2j]))-ylin)**2).sum()) for w1i in w1range] for w2j in w2range])
to:
cost = np.array([[np.log(((X.dot(np.array([w1i,w2j]))-y)**2).sum()) for w1i in w1range] for w2j in w2range])
Changed lines 125-133 from:
to:
De manière optionnelle, on peut essayer de tracer la courbe de coût en 3D pour mieux visualiser ce qui se passe:

(:source lang=python:)
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(w1, w2, cost, rstride = 1, cstride=1 )
(:sourceend:)
Changed line 39 from:
* {$\hat b = E(Y) - cov(X, Y)/\sigma^2_x E(X)$}
to:
* {$\hat b = E(Y) - \frac{cov(X, Y)}{\sigma^2_x} E(X)$}
Changed line 36 from:
!!!! Estimation de paramètres probabiliste
to:
!!!! Estimation de paramètres probabilistes
Changed line 14 from:
* {$y = ax + b + \epsilon,\ \epsilon \sim \mathcal N(0,\sigma)$}
to:
* {$y = ax + b + \epsilon,\ \epsilon \sim \mathcal N(0,\sigma^2)$}
Changed line 22 from:
sig = .4
to:
sig = .4 # écart type
Added lines 171-173:
* Construire plusieurs métriques d'évaluation
** Erreur moyenne de prédiction
** En arrondissant les sorties de notre modèle, on peut compter le nombre de fois où l'on trouve la note exacte
Changed lines 149-150 from:
* Format pkl: [[(Attach:)wine.pkl]]
to:

!!! Chargement des données

Le jeu de données sera séparé en un jeu d'apprentissage
(pour l'estimation des paramètres) et un jeu de test (pour l'estimation non biaisée de la performance)

(:source lang=python:)
data = np
.loadtxt("ressources/winequality-red.csv", delimiter=";", skiprows=1)
N,d = data.shape # extraction des dimensions
pcTrain  = 0.7 # 70% des données en apprentissage
allindex = np.random.permutation(N)
indTrain = allindex[:int(pcTrain*N)
]
indTest = allindex[int(pcTrain*N):]
X = data[indTrain,:-1] # pas la dernière colonne (= note à prédire)
Y = data[indTrain,-1]  # dernière colonne (= note à prédire)
# Echantillon de test (pour la validation des résultats)
XT = data[indTest,:-1] # pas la dernière colonne (= note à prédire)
YT = data[indTest,-1]  # dernière colonne (= note à prédire)
(:sourceend:)

Added lines 169-172:
* Utiliser le modèle analytique (qui va beaucoup plus vite sur les petits jeux de données)
* Toujours comparer les performances en apprentissage et en test
* Construire de nouveau attribut (ie de nouvelles colonnes descriptives) en combinant les données existantes (à la manière de ce que nous avons fait précédemment)

Changed lines 68-86 from:
!!!%subsection% Extension non-linéaire

Générer de nouvelles données non-linéaires selon la formule: {$yquad = ax^2 + b x + c + \epsilon,\ \epsilon \sim \mathcal N(0,\sigma)$} avec les mêmes valeurs de paramètres que précédemment.

%width=300px% Attach:TME10_datagenNL.png

Estimation de {$a,b,c$} au sens de l'erreur des moindres carrés:
* construire un nouveau {$Xe= [x^2,\ x, 1]$} en utilisant la méthode @@hstack@@ vue dans la question précédente
* estimer un vecteur w sur ces données
* reconstruire {$yquad$}
* calculer l'erreur moyenne de reconstruction


Vous devez obtenir un estimateur de la forme:

%width=300px% Attach:TME10_quad.png


to:
Changed lines 125-167 from:
!!%section% Données réelles
to:

!!!%subsection% Extension non-linéaire

Générer de nouvelles données non-linéaires selon la formule: {$yquad = ax^2 + b x + c + \epsilon,\ \epsilon \sim \mathcal N(0,\sigma)$} avec les mêmes valeurs de paramètres que précédemment.

%width=300px% Attach:TME10_datagenNL.png

Estimation de {$a,b,c$} au sens de l'erreur des moindres carrés:
* construire un nouveau {$Xe= [x^2,\ x, 1]$} en utilisant la méthode @@hstack@@ vue dans la question précédente
* estimer un vecteur w sur ces données
* reconstruire {$yquad$}
* calculer l'erreur moyenne de reconstruction


Vous devez obtenir un estimateur de la forme:

%width=300px% Attach:TME10_quad.png


!!%section% Données réelles

Récupérer les données sur la qualité vin.

* Source originale: [[https://archive.ics.uci.edu/ml/datasets/Wine+Quality]]
* Format pkl: [[(Attach:)wine.pkl]]

Tester différents modèles de régression sur ce jeu de données.

!!%section% Pour aller plus loin

!!!%subsection% Optimisation par descente de gradient '''stochastique'''

Soit un problème avec des données {$(x_i,y_i)_{i=1,\ldots,N}$}, une fonction de décision/prédiction paramétrée par un vecteur {$w$} et une fonction de cout à optimiser {$C(w)$}.


l'algorithme de la descente de gradient '''stochastique''' est le suivant :

# {$w_0 \leftarrow init$} par exemple : 0
# boucle
## tirage d'une donnée {$i$}: {$(x_i,y_i)$}
## {$w_{t+1} \leftarrow w_{t} - \epsilon \nabla_w C_i(w)$}

Reprendre rapidement les questions ci-dessus.
Changed lines 99-152 from:
Plusieurs critères d'arrêt sont envisageables: nombre d'itérations, plus de mouvement sur {$w$}, norme de {$\nabla_w C(w)$} suffisamment petite...

'''Q1''' : programmer cette méthode et comparer les résultats avec la solution analytique

Visualiser la solution obtenue:
  % soit x et y et une solution (a,b)
  figure
  plot(x,y,'b+');  % affichage des points
  hold on; 
       % ne pas écraser l'affichage précédent
  yhat = a*x+b;
  plot(x,yhat,'r'); % sortie du modèle

NB
: il est possible de visualiser l'évolution du modèle au cours des itérations du gradient

  % dans la boucle du
gradient, en notant (at,bt) la solution courante et avec une figure déjà ouverte
  hold off % pour que l'affichage écrase la figure précédente
  plot
(x,y,'b+');  % affichage des points
  hold on;          % ne pas écraser l
'affichage précédent
 
plot(x,at*x+bt,'r');
  drawnow           % forcer l'affichage maintenant
  pause(1);        % 1 seconde de pause pour laisser le temps de voir qqch

'''Q2'''
: Dans le cas de la régression 1D, il n'y a que deux coefficients {$(a,b) = (w_1,w_2)$}. Dans l'espace de ces paramètres, tracer l'évolution de la solution au fil des itérations. Vous placerez le point correspondant à la solution analytique pour analyser la qualité de la solution obtenue

'''Q3''' : Etudier la qualité de la solution en fonction de {$\epsilon$}. Vous sauverez 3 images correspondant aux solutions {$\epsilon$} bon, trop petit ou trop grand.

'''Q4''' : Afin d'encore mieux comprendre ce qui se passe, je vous propose d'ajouter un élément d'affichage: la fonction coût.
Pour l'instant, la solution analytique nous donne le minimum mais on peut en fait dessiner cette fonction. Le processus (non trivial) est décrit ci-dessous:

# Génération d'un ensemble de couple {$(a_j,b_j)_{j=1,\ldots,G}$} maillant l'ensemble de l'espace
# Evaluation du coût pour tous ces couples
# Extraction et affichage des isocontours de coût pour visualiser C

Eléments de code pour l'implémentation:

  % génération des couples
  %- plage des pour a
  a = -2:0.1:5;
  %- plage des pour b
  b = -4:0.1:10;
  %- génération des couples: fonction meshgrid
  [w1,w2]=meshgrid(a,b); % toutes les combinaisons sont générées
 
  % Evaluation du cout pour tous les couples
  ... % A FAIRE après avoir noté la forme des matrices w1,w2
      % A stocker dans une matrice C de même dimension que w1 et w2

  % affichage
  contour(w1,w2,C);

On cherche un résultat de la forme:
%center% %width=300px% Attach:plotCost.png

to:
Compléter le squelette d'implémentation fourni ci-dessous:

(:source lang=python:)
wstar = np.linalg.solve(X.T.dot(X), X.T.dot(ylin)) # pour se rappeler du
w optimal

eps = 5e-3
nIterations = 30

w = np.zeros(X.shape[1]) # init à 0
allw = [w]
for i in xrange(nIterations):
    # A COMPLETER => calcul du gradient vu en TD
    allw.append
(w)
    print w

allw = np.array
(allw)
(:sourceend:)

* Tester différentes valeurs d'epsilon
* Tester différentes initialisations
* comparer les résultats théoriques (solution analytique) et par descente de gradient

On s'intéresse ensuite à comprendre la descente de gradient dans l'espace des paramètres. Le code ci-dessous permet de
:
* Tracer le cout pour un ensemble de paramètres
* Comprendre comment la solution évolue (en terme de coût) au fil des itérations du
gradient

(:source lang=python:)
# tracer de l'espace des couts
ngrid = 20
w1range = np.linspace(-0.5, 8, ngrid)
w2range = np.linspace(-1.5, 1.5, ngrid)
w1,w2 = np.meshgrid(w1range,w2range)

cost = np.array
([[np.log(((X.dot(np.array([w1i,w2j]))-ylin)**2).sum()) for w1i in w1range] for w2j in w2range])

plt.figure()
plt.contour(w1, w2, cost)
plt.scatter(wstar[0], wstar[1],c=
'r')
plt.
plot(allw[:,0],allw[:,1],'b+-' ,lw=2 )
(:sourceend:)

Vous devez obtenir une figure :

%width=300px% Attach
:TME10_gradient.png
Added lines 84-152:



!!!%subsection% Optimisation par descente de gradient

Soit un problème avec des données {$(x_i,y_i)_{i=1,\ldots,N}$}, une fonction de décision/prédiction paramétrée par un vecteur {$w$} et une fonction de cout à optimiser {$C(w)$}.
Notre but est de trouver les paramètres {$w^\star$} minimisant la fonction de coût:
{$$ w^\star = \arg\min_w C(w)$$}

l'algorithme de la descente de gradient est le suivant (rappel):

# {$w_0 \leftarrow init$} par exemple : 0
# boucle
## {$w_{t+1} \leftarrow w_{t} - \epsilon \nabla_w C(w)$}

Plusieurs critères d'arrêt sont envisageables: nombre d'itérations, plus de mouvement sur {$w$}, norme de {$\nabla_w C(w)$} suffisamment petite...

'''Q1''' : programmer cette méthode et comparer les résultats avec la solution analytique

Visualiser la solution obtenue:
  % soit x et y et une solution (a,b)
  figure
  plot(x,y,'b+');  % affichage des points
  hold on;          % ne pas écraser l'affichage précédent
  yhat = a*x+b;
  plot(x,yhat,'r'); % sortie du modèle

NB: il est possible de visualiser l'évolution du modèle au cours des itérations du gradient

  % dans la boucle du gradient, en notant (at,bt) la solution courante et avec une figure déjà ouverte
  hold off % pour que l'affichage écrase la figure précédente
  plot(x,y,'b+');  % affichage des points
  hold on;          % ne pas écraser l'affichage précédent
  plot(x,at*x+bt,'r');
  drawnow          % forcer l'affichage maintenant
  pause(1);        % 1 seconde de pause pour laisser le temps de voir qqch

'''Q2''' : Dans le cas de la régression 1D, il n'y a que deux coefficients {$(a,b) = (w_1,w_2)$}. Dans l'espace de ces paramètres, tracer l'évolution de la solution au fil des itérations. Vous placerez le point correspondant à la solution analytique pour analyser la qualité de la solution obtenue

'''Q3''' : Etudier la qualité de la solution en fonction de {$\epsilon$}. Vous sauverez 3 images correspondant aux solutions {$\epsilon$} bon, trop petit ou trop grand.

'''Q4''' : Afin d'encore mieux comprendre ce qui se passe, je vous propose d'ajouter un élément d'affichage: la fonction coût.
Pour l'instant, la solution analytique nous donne le minimum mais on peut en fait dessiner cette fonction. Le processus (non trivial) est décrit ci-dessous:

# Génération d'un ensemble de couple {$(a_j,b_j)_{j=1,\ldots,G}$} maillant l'ensemble de l'espace
# Evaluation du coût pour tous ces couples
# Extraction et affichage des isocontours de coût pour visualiser C

Eléments de code pour l'implémentation:

  % génération des couples
  %- plage des pour a
  a = -2:0.1:5;
  %- plage des pour b
  b = -4:0.1:10;
  %- génération des couples: fonction meshgrid
  [w1,w2]=meshgrid(a,b); % toutes les combinaisons sont générées
 
  % Evaluation du cout pour tous les couples
  ... % A FAIRE après avoir noté la forme des matrices w1,w2
      % A stocker dans une matrice C de même dimension que w1 et w2

  % affichage
  contour(w1,w2,C);

On cherche un résultat de la forme:
%center% %width=300px% Attach:plotCost.png

Changed lines 32-33 from:
!!%section% Validation des formules analytiques
to:
!!!%subsection% Validation des formules analytiques
Changed lines 36-37 from:
!!! Estimation de paramètres probabiliste
to:
!!!! Estimation de paramètres probabiliste
Changed lines 45-46 from:
!!! Formulation au sens des moindres carrés
to:
!!!! Formulation au sens des moindres carrés
Changed lines 68-69 from:
!!%section% Extension non-linéaire
to:
!!!%subsection% Extension non-linéaire
Added lines 79-85:


Vous devez obtenir un estimateur de la forme:

%width=300px% Attach:TME10_quad.png

!!%section% Données réelles
Changed line 72 from:
Attach:TME10_datagenNL.png
to:
%width=300px% Attach:TME10_datagenNL.png
Changed lines 70-71 from:
Générer de nouvelles données non-linéaires selon la formule: {$y = ax^2 + b^x + c + \epsilon,\ \epsilon \sim \mathcal N(0,\sigma)$} avec les mêmes valeurs de paramètres que précédemment.
to:
Générer de nouvelles données non-linéaires selon la formule: {$yquad = ax^2 + b x + c + \epsilon,\ \epsilon \sim \mathcal N(0,\sigma)$} avec les mêmes valeurs de paramètres que précédemment.
Added lines 74-78:
Estimation de {$a,b,c$} au sens de l'erreur des moindres carrés:
* construire un nouveau {$Xe= [x^2,\ x, 1]$} en utilisant la méthode @@hstack@@ vue dans la question précédente
* estimer un vecteur w sur ces données
* reconstruire {$yquad$}
* calculer l'erreur moyenne de reconstruction
Changed lines 66-73 from:
Que se passe-t-il si on n'ajoute pas de biais dans {$X$}? Tracer la solution obtenue.
to:
Que se passe-t-il si on n'ajoute pas de biais dans {$X$}? Tracer la solution obtenue.

!!%section% Extension non-linéaire

Générer de nouvelles données non-linéaires selon la formule: {$y = ax^2 + b^x + c + \epsilon,\ \epsilon \sim \mathcal N(0,\sigma)$} avec les mêmes valeurs de paramètres que précédemment.

Attach:TME10_datagenNL.png

Changed line 22 from:
sig = .3
to:
sig = .4
Changed lines 43-66 from:
%width=300px% Attach:TME10_linEstP.png
to:
%width=300px% Attach:TME10_linEstP.png

!!! Formulation au sens des moindres carrés

Nous partons directement sur une écriture matricielle. Du coup, il est nécessaire de construire la matrice {$X$}... Il y a plein de piège (mais je donne la solution plus bas):
* usage de @@hstack@@ pour mettre cote à cote deux matrices: ATTENTION, la fonction prend en argument un tuple
* obligation de faire un reshape sur x (il ne sait pas comment agréger une matrice et un vecteur)

Solution:
(:source lang=python:)
X = np.hstack((x.reshape(N,1),np.ones((N,1))))
(:sourceend:)

Il faut ensuite poser et résoudre un système d'équations linéaires de la forme {$A w = B$}

Rappel des formules vues en cours/TD:
* {$A = X^T X$} (produit matriciel)
*{$B = X^T Y$} (produit matriciel)

Fonction de résolution: @@np.linalg.solve(A,B)@@


Vous devez obtenir la même solution que précédemment.
Que se passe-t-il si on n'ajoute pas de biais dans {$X$}? Tracer la solution obtenue.
Changed lines 38-42 from:
* {$\hat a=\cov(X,Y) / \sigma_x^2$}
* {$\hat b = E(Y) - \cov(X, Y)/\sigma^2_x E(X)$}

Estimer les paramètres
to:
* {$\hat a = cov(X,Y) / \sigma_x^2$}
* {$\hat b = E(Y) - cov(X, Y)/\sigma^2_x E(X)$}

Estimer les paramètres, faire passer les {$x$} dans le modèle et tracer la courbe associée.

%width=300px% Attach:TME10_linEstP.png
Changed line 13 from:
* {$x\in [0,\ 1]$} (tirage avec un simple @@rand()@@
to:
* {$x\in [0,\ 1]$} tirage avec un simple @@rand()@@
Added lines 15-24:
* Pour faire cet énoncé, je suis parti des valeurs suivantes:

(:source lang=python:)
a = 6.
b = -1.
c = 1
N = 100
sig = .3
(:sourceend:)

Added lines 31-41:

!!%section% Validation des formules analytiques

Nous avons vu deux types de résolutions analytique: à partir des estimateurs des espérances et co-variances d'une part et des moindres carrés d'autre part. Testons les deux méthodes.

!!! Estimation de paramètres probabiliste

* {$\hat a=\cov(X,Y) / \sigma_x^2$}
* {$\hat b = E(Y) - \cov(X, Y)/\sigma^2_x E(X)$}

Estimer les paramètres
Changed line 13 from:
* {$x\in [0,\1]$} (tirage avec un simple @@rand()@@
to:
* {$x\in [0,\ 1]$} (tirage avec un simple @@rand()@@
Changed lines 11-21 from:
Dans un premier temps, générons des données jouets
to:
Dans un premier temps, générons des données jouets paramétriques:
* {$N$} nombre de points à générer
* {$x\in [0,\1]$} (tirage avec un simple @@rand()@@
* {$y = ax + b + \epsilon,\ \epsilon \sim \mathcal N(0,\sigma)$}
Notes:
* Vous utiliserez la fonction @@np.random.randn@@ pour le tirage de données selon une loi normale
* Vous obtenez un écart-type {$\sigma$} pas simple multiplication sur le tirage précédent
** il s'agit en fait d'une homothétie... Si 66% des points étaient entre -1 et 1, en multipliant par {$\sigma$}, 66% des points sont entre {$-\sigma$} et {$\sigma$}

%width=300px% Attach:TME10_datagen.png

Added lines 1-11:
%define=subsection color=purple%
%define=section color=blue%


! TME 10: Régression

Nous travaillons aujourd'hui sur les méthodes de régression. Dans un premier temps, nous allons générer des données jouets pour valider l'implémentation des formules vues en cours. Nous verrons ensuite comment passer au non-linéaire (toujours sur des données jouets) puis nous étudierons ces méthodes sur des données réelles (qualité du vin).

!!%section% Données jouets

Dans un premier temps, générons des données jouets