Cours

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

Données jouets

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^2)$}
  • Pour faire cet énoncé, je suis parti des valeurs suivantes:
a = 6.
b = -1.
N = 100
sig = .4 # écart type

Notes:

  • Vous utiliserez la fonction np.random.randn pour le tirage de données selon une loi normale
  • Vous obtenez un écart-type {$\sigma$} par 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$}

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 probabilistes

  • {$\hat a = cov(X,Y) / \sigma_x^2$}
  • {$\hat b = E(Y) - \frac{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.

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:

X = np.hstack((x.reshape(N,1),np.ones((N,1))))

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.

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):

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

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

wstar = np.linalg.solve(X.T.dot(X), X.T.dot(y)) # 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)
  • 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
# 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]))-y)**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 )

Vous devez obtenir une figure :

De manière optionnelle, on peut essayer de tracer la courbe de coût en 3D pour mieux visualiser ce qui se passe:

from mpl_toolkits.mplot3d import Axes3D

costPath = np.array([np.log(((X.dot(wtmp)-y)**2).sum()) for wtmp in allw])
costOpt  = np.log(((X.dot(wstar)-y)**2).sum())

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(w1, w2, cost, rstride = 1, cstride=1 )
ax.scatter(wstar[0], wstar[1],costOpt, c='r')
ax.plot(allw[:,0],allw[:,1],costPath, 'b+-' ,lw=2 )

Extension non-linéaire (solution analytique)

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.

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 en utilisant la solution analytique (les données sont petites, pas la peine de s'embéter)
  • reconstruire {$yquad$}
  • calculer l'erreur moyenne de reconstruction

Vous devez obtenir un estimateur de la forme:

Données réelles

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

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)

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)

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

  • 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 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
  • 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)

Pour aller plus loin

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 :

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

Reprendre rapidement les questions ci-dessus.