Evaluation numérique du risque: Illustrations

Le but de ce TP est de voir comment, sur un exemple simple, on peut evaluer en pratique le risque moyen ou de généralisation.

Nous allons nous mettre dans une situation où on est capable de tout simuler facilement, c'est à dire qu'on va se fixer la loi jointes des entrées/sorties $(X,Y)$, et on pourra donc

• évaluer le risque de Bayes
• estimer l'erreur de généralisation
• estimer l'erreur moyenne de généralisation.

Dans ce TP on va utiliser la librairie scikit-learn.

https://scikit-learn.org/stable/index.html

Définition de la loi jointe et estimation du risque de Bayes

Introduisons maintenant la situation que nous allons considérer. On va considérer une regression non linéaire simple.

$$ Y = \arctan(\beta_0+\beta_1 X) +\varepsilon$$

où $\varepsilon\sim N(0,0.2^2)$ et $X\sim N(0,10^2)$. Ces relations définissent la loi jointe de $(X,Y)$.

import matplotlib.pyplot as plt
import numpy as np
import numpy.random as rnd

Fixons nous $(\beta_0,\beta_1)$. On va les choisir aléatoirement, mais ceci n'a rien à voir avec la loi jointes des entrées/sorties.

p = 1
beta = rnd.normal(0, 1, p + 1)
sig = 0.2

beta

array([-0.63770377, -1.0142289 ])

On va se donner $n$ observations

n = 50
x = rnd.normal(0, 10, (n, p))

Vous pouvez rajouter des dimensions aux entrées, mais ces moins parlant pour des illustrations graphiques.

X = np.concatenate((np.ones((n,1)), x), axis=1)
X.shape

(50, 2)

Y = np.arctan(X.dot(beta)) + rnd.normal(0, sig, n)

On peut alors estimer le risque de Bayes avec comme fonction de coût le coût quadratique:

$$L(y,y') = (y-y')^2.$$

Comme on l'a vu en cours, le risque de Bayes est donné ici par $R^\ast=Var(\varepsilon)=\sigma^2 = 0.2^2=0.04$ et le predicteur de Bayes est bien entendu

$$\eta(X) = \arctan(\beta_0+\beta_1 X).$$

Maintenant illustrons nos données avec le prédicteur de Bayes.

n0 = 100
x_plot = np.linspace(x.min(), x.max(), n0)
pred_bayes = np.arctan(beta[0] + beta[1]*x_plot)

plt.plot(x, Y, 'o', label='data')
plt.plot(x_plot, pred_bayes, label='predict Bayes')
plt.legend()
plt.show()

Définition de la classe de modèle utilisée pour la regression

On est dans un context non linéaire, et on va considérer une regression à noyau. Vous pourrait tester d'autres noyaux et/ou modèles (k-NN, forêts aléatoire, réseaux de neuronne...), mais ici on se contentera d'un noyau Gaussien dont on sait qu'il est universel. On importe ce type de modèle depuis sci-kit learn

from sklearn.kernel_ridge import KernelRidge

On fixe "à la louche" un paramètre de régularisation alpha et la variance du noyau Gaussien gamma. On verra plus bas comment les choisir en pratique.

model = KernelRidge(alpha=0.06, kernel='rbf', gamma=0.01)

Minimiseur du risque empirique.

On cherche une solution au problème de minimisation sous contrainte avec le risque empirique avec le noyau considéré (on avait vu en cours comment déterminer explicitement la solution dans ce cas). Avec scikit-learn il suffit d'utiliser la méthode fit.

model.fit(x, Y)

KernelRidge(alpha=0.06, gamma=0.01, kernel='rbf')

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

On a donc obtenu $\hat f _S$, qui dépend bien entendu du jeu d'exemple généré plus haut. Dans notre cas très simple on peut tracer $\hat f_S$, en utilisant la méthode predict pour calculer les images des points (qui sont en effet des predictions pour les points considérés).

pred = model.predict(x_plot.reshape((n0,1)))

plt.plot(x, Y, 'o', label='data')
plt.plot(x_plot, pred_bayes, label='pred Bayes')
plt.plot(x_plot, pred, label='predictions')
plt.legend()
plt.show()

Risque de généralisation et risque moyen

Estimons maintenant l'erreur/risque de généralisation de $\hat f_S$, qui je le rappelle dépend des exemples générés plus haut, et l'erreur/risque moyen. Nous allons procéder à l'aide de la loi des grands nombres. Comme on considère un coû quadratique, les risques empiriques pourront être calculé à l'aide de la MSE (mean-square-error) qu'on peut importer de scikit-learn.

from sklearn.metrics import mean_squared_error as MSE

Pour le risque de généralisation, on génère de nouvelles donées.

n_estim = 20000
x_new = rnd.normal(0, 10, (n_estim, p))
X_new = np.concatenate((np.ones((n_estim,1)), x_new), axis=1)
Y_new = np.arctan(X_new.dot(beta)) + rnd.normal(0, sig, n_estim)
Y_pred = model.predict(x_new.reshape((n_estim,1)))

et on a pour risque de généralisation

MSE(Y_new,Y_pred)

0.10814719322695762

Pour obtenir le risque moyen, on doit prendre l'espérance du risque de généralisation sur les exemples. On va donc devoir faire une boucle pour générer des exemples, calculer le predicteur qui minimise le risque empirique, et évaluer l'erreur de généralisation. On va prendre des valeurs assez faible pour faire les estimations des espérances pour ne pas faire tourner la machine trop longtemps.

m_estim = 20
n_estim = 2000
mse = 0.

for _ in range(m_estim):
    x_new = rnd.normal(0, 10, (n_estim, p))
    X_new = np.concatenate((np.ones((n_estim,1)), x_new), axis=1)
    Y_new = np.arctan(X_new.dot(beta)) + rnd.normal(0, sig, n_estim)
    model = KernelRidge(alpha=0.06, kernel='rbf', gamma=0.01)
    model.fit(x_new, Y_new)
    Y_pred = model.predict(x_new.reshape((n_estim,1)))
    mse += MSE(Y_new,Y_pred)

le risque moyen est donc donnée par

mse/m_estim

0.05581331855816599

On voit que l'estimation du risque moyen est plutôt bonne et que le risque de généralisation est plus élevé. En effet, pour ce dernier, le prédicteur dépend du jeu d'exemple.

Effect du paramètre de régularisation

On va s'intéresser à l'erreur de généralisation en fonction du paramètre de régularisation alpha. On avait vu en cours l'effet potentiel de ce paramètre sur la qualité de prédiction.

n_estim = 20000
risque_alpha = []
range_alpha = np.arange(0.0001, 0.2, 0.005)

for alpha in range_alpha:
    model = KernelRidge(alpha=alpha, kernel='rbf', gamma=0.01)
    model.fit(x, Y)
    x_new = rnd.normal(0, 10, (n_estim, p))
    X_new = np.concatenate((np.ones((n_estim,1)), x_new), axis=1)
    Y_new = np.arctan(X_new.dot(beta)) + rnd.normal(0, sig, n_estim)
    Y_pred = model.predict(x_new.reshape((n_estim,1)))
    risque_alpha.append(MSE(Y_new,Y_pred))

alpha_opt = range_alpha[np.argmin(risque_alpha)]
alpha_opt

0.0051

risque_alpha[np.argmin(risque_alpha)]

0.0963244389375491

Voici l'évolution du risque de généralisation en fonction du paramètre de régularisation.

plt.plot(range_alpha, risque_alpha)
plt.show()

On retrouve l'idée de la courbe rouge de la figure 7.1 p.220 de la deuxième référence. Si on change les données d'exemple cette courbe sera différente! Même si ici nous n'avons pas beaucoup de paramètre, on voit le rôle du paramètre de régularisation. Plus alpha est petit et plus notre modèle peut être compliqué, sur-apprendre, et donc mal prédire. D'où un risque de généralisation élevé. Aussi, plus alpha augmente, plus la compléxité du modèle diminue, le modèle devient trop simple, il va sous apprendre, et même mal prédire. C'est ce qu'on peut voir sur les deux figures suivantes.

model = KernelRidge(alpha=0.0000001, kernel='rbf', gamma=0.01)
model.fit(x, Y)

x_plot = np.linspace(x.min(), x.max(), 100)
pred = model.predict(x_plot.reshape((-1,1)))
pred_bayes = np.arctan(beta[0] + beta[1]*x_plot)

plt.plot(x, Y, 'o', label='data', alpha=0.2)
plt.plot(x_plot, pred_bayes, label='pred Bayes')
plt.plot(x_plot, pred, label='predictions')
plt.legend()
plt.show()

Le modèle cherche vraiment à sur-apprendre les données d'entrainements.

model = KernelRidge(alpha=2, kernel='rbf', gamma=0.01)
model.fit(x, Y)

x_plot = np.linspace(x.min(), x.max(), 100)
pred = model.predict(x_plot.reshape((-1,1)))
pred_bayes = np.arctan(beta[0] + beta[1]*x_plot)

plt.plot(x, Y, 'o', label='data', alpha=0.2)
plt.plot(x_plot, pred_bayes, label='pred Bayes')
plt.plot(x_plot, pred, label='predictions')
plt.legend()
plt.show()

Ici le modèle sous-apprend et a même du mal à coller aux données d'entrainements. Enfin, en prenant un alpha qui minimise l'erreur de généralisation on a un modèle de compléxité intermédiaire. Mais nous verrons dans la section suivante comment fixer les paramètres de notre modèle alpha et gamma, appelés hyperparamètres.

model = KernelRidge(alpha=alpha_opt, kernel='rbf', gamma=0.01)
model.fit(x, Y)

x_plot = np.linspace(x.min(), x.max(), 100)
pred = model.predict(x_plot.reshape((-1,1)))
pred_bayes = np.arctan(beta[0] + beta[1]*x_plot)

plt.plot(x, Y, 'o', label='data', alpha=0.2)
plt.plot(x_plot, pred_bayes, label='pred Bayes')
plt.plot(x_plot, pred, label='predictions')
plt.legend()
plt.show()

Cross validation

Comment fait-on en pratique pour évaluer ces risques sachant

• qu'on ne connait pas la loi jointe, et qu'on ne peut pas les estimer comme on vient de le faire,
• qu'on a qu'un jeu d'exemple fixe à disposition.

On remarquera qu'on est plus intéressé par le risque de généralisation, puisque c'est lui qui mesure vraiment les erreurs de généralisation qu'on fera en pratique.

Une technique pour se faire est la cross validation. Je vous renvoie aux références 2 et 3 que je vous avez donné en début de cours pour les subtilités de la technique, et à la page de scikit-learn qui possède cette fonctionnalité. Il y a aussi une page wikipedia que vous pouvez consulter.

https://scikit-learn.org/stable/modules/cross_validation.html

On va se donner une jeu d'exemple plus gros que précédement, pour avoir un découpage avec des données de tests par trop petit pour la cross validation.

n = 150
x = rnd.normal(0, 10, (n, p))
X = np.concatenate((np.ones((n,1)), x), axis=1)
Y = np.arctan(X.dot(beta)) + rnd.normal(0, sig, n)

model = KernelRidge(alpha=0.06, kernel='rbf', gamma=0.01)
model.fit(x, Y)

x_plot = np.linspace(x.min(), x.max(), 100)
pred = model.predict(x_plot.reshape((-1,1)))
pred_bayes = np.arctan(beta[0] + beta[1]*x_plot)

plt.plot(x, Y, 'o', label='data', alpha=0.2)
plt.plot(x_plot, pred_bayes, label='pred Bayes')
plt.plot(x_plot, pred, label='predictions')
plt.legend()
plt.show()

Avec ce nouveau jeu de données, on peut estimer le risque de généralisation

n_estim = 20000
x_new = rnd.normal(0, 10, (n_estim, p))
X_new = np.concatenate((np.ones((n_estim,1)), x_new), axis=1)
Y_new = np.arctan(X_new.dot(beta)) + rnd.normal(0, sig, n_estim)
Y_pred = model.predict(x_new.reshape((n_estim,1)))

MSE(Y_new,Y_pred)

0.07029323368331582

On peut effectuer la cross validation (de base) très facilement avec sci-kit learn.

from sklearn.model_selection import cross_val_score

model = KernelRidge(alpha=0.06, kernel='rbf', gamma=0.01)
scores = cross_val_score(model, x, Y, cv=5, scoring='neg_mean_squared_error')

-scores.mean()

0.08580735407908234

Dans notre cas très simple l'approximation, par cross validation, des risques moyen et de généralisation est plutôt bonne. Je rappelle que le noyau Gaussien est universel!

Recherche des hypers paramètres

Il s'agit de trouver les paramètres alpha et gamma qui vont minimiser le risque. Quel risque? Celui qu'on peut estimer, par cross validation par exemple. Cela nous permet d'éviter de sur- ou sous-apprendre. Pour cela on peut utiliser la fonction GridSearchCV.

from sklearn.model_selection import GridSearchCV

parameters = {'alpha': np.arange(0.0001, 0.2, 0.005),
              'gamma': np.arange(0.0001, 0.2, 0.005)}
kr = KernelRidge(kernel='rbf')      
model = GridSearchCV(kr, parameters, cv=3, scoring='neg_mean_squared_error')

model.fit(x, Y)

GridSearchCV(cv=3, estimator=KernelRidge(kernel='rbf'),
             param_grid={'alpha': array([1.000e-04, 5.100e-03, 1.010e-02, 1.510e-02, 2.010e-02, 2.510e-02,
       3.010e-02, 3.510e-02, 4.010e-02, 4.510e-02, 5.010e-02, 5.510e-02,
       6.010e-02, 6.510e-02, 7.010e-02, 7.510e-02, 8.010e-02, 8.510e-02,
       9.010e-02, 9.510e-02, 1.001e-01, 1.051e-01, 1.101e-01, 1.151e-01,
       1.201e-01, 1.251e-01, 1.301e-01...
       3.010e-02, 3.510e-02, 4.010e-02, 4.510e-02, 5.010e-02, 5.510e-02,
       6.010e-02, 6.510e-02, 7.010e-02, 7.510e-02, 8.010e-02, 8.510e-02,
       9.010e-02, 9.510e-02, 1.001e-01, 1.051e-01, 1.101e-01, 1.151e-01,
       1.201e-01, 1.251e-01, 1.301e-01, 1.351e-01, 1.401e-01, 1.451e-01,
       1.501e-01, 1.551e-01, 1.601e-01, 1.651e-01, 1.701e-01, 1.751e-01,
       1.801e-01, 1.851e-01, 1.901e-01, 1.951e-01])},
             scoring='neg_mean_squared_error')

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Les meilleurs paramètres trouvés sont

model.best_params_

{'alpha': 0.10010000000000001, 'gamma': 0.0101}

x_plot = np.linspace(x.min(), x.max(), 100)
pred = model.predict(x_plot.reshape((-1,1)))
pred_bayes = np.arctan(beta[0] + beta[1]*x_plot)

plt.plot(x, Y, 'o', label='data', alpha=0.2)
plt.plot(x_plot, pred_bayes, label='pred Bayes')
plt.plot(x_plot, pred, label='predictions')
plt.legend()
plt.show()

scikit-learn propose tout un tas de fonctionnalité (optimisée) permettant d'automatiser l'évaluation des modèles sur des jeu de donnée, mais je vous laisse vous reporter à la documentation.

Enfin, on peut regarder le modèle d'erreur. Ici, c'est très simple. Si on retranche les predictions aux points x, on doit retrouver à peu près une erreur Gaussiennne centrée de variance $0.2^2$.

err = Y-model.predict(x)

Pour tester la normalité de err on va utiliser un test statistique, et si on peut considérer comme un échantillon Gaussien on estimera ces paramètres.

from scipy.stats import normaltest
from scipy.stats import norm

normaltest(err)

NormaltestResult(statistic=2.5079239158887283, pvalue=0.28537192258619043)

La $p$-value est grande, on ne peut donc pas rejeter le fait qu'il s'agit d'un échantillon Gaussien. Estimons ces paramètres.

norm.fit(err)

(-0.00042547455007194184, 0.2709213153648689)

On retrouve une erreur centrée, avec un ecart type proche du vraie sig=0.2. On peut enfin construire un intervalle de confiance pour les prédictions.

m_err, sig_err = norm.fit(err)

q = norm.ppf(1-0.05/2) # quantile de la loi N(0,1).

plt.plot(x, Y, 'o', label='data', alpha=0.2)
plt.plot(x_plot, pred_bayes, label='pred Bayes')
plt.plot(x_plot, pred, label='predictions')
plt.fill_between(x_plot, pred-sig_err*q, pred+sig_err*q, color="g", alpha=0.1, label='IC 95%')
plt.legend()
plt.show()