{ "metadata": { "kernelspec": { "name": "python", "display_name": "Python (Pyodide)", "language": "python" }, "language_info": { "codemirror_mode": { "name": "python", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8" } }, "nbformat_minor": 4, "nbformat": 4, "cells": [ { "cell_type": "markdown", "source": "# TP 1 – Pour aller plus loin", "metadata": {} }, { "cell_type": "markdown", "source": "### 1 : Régression linéaire univariée sans scikit-learn", "metadata": {} }, { "cell_type": "markdown", "source": "Notre but ici va être d'inférer le poids d'un manchot à partir de la longueur de sa palette natatoire. Nous sommes donc dans le cas particulier où $p=1$.", "metadata": {} }, { "cell_type": "code", "source": "(X_train, X_test, y_train, y_test) = model_selection.train_test_split(X, y_regress, \n test_size=0.3, random_state=25)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "x_3 = X_train[:, -1] # on ne garde que la dernière colonne/variable de X", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "plt.scatter(x_3, y_train)\nplt.xlabel(\"Longueur de la palette natatoire (mm)\")\nplt.ylabel(\"Poids (g)\")", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "Il s'agit donc d'apprendre une fonction $f: x \\mapsto a x + b$ telle que $y \\approx f(x_3)$.\n\nCe problème se résout grâce à la méthode des moindres carrés : il s'agit de trouver $a, b$ qui minimisent $J(a, b) = \\frac{1}{n} \\sum_{i=1}^n ((a x_i + b) - y_i)^2$. \n\nIl s'agit d'une forme quadratique, convexe, dont on trouve le minimum en annulant sa dérivée en $a$ et en $b$.\n\n$\\frac{\\partial J}{\\partial a} (a, b)= \\frac2n \\sum_{i=1}^n x_i^2 a + \\frac2n \\sum_{i=1}^n x_i b - \\frac2n \\sum_{i=1}^n x_i y_i$\n\net\n\n$\\frac{\\partial J}{\\partial b} (a, b)= b + \\frac2n \\sum_{i=1}^n x_i a - \\frac2n \\sum_{i=1}^n y_i$\n\nCe dont on déduit, à condition que $\\frac1n \\sum_{i=1}^n x_i^2 \\neq \\left(\\frac1n \\sum_{i=1}^n x_i\\right)^2$,\n\n$a = \\frac{\\frac1n \\sum_{i=1}^n x_i y_i - \\frac1n \\sum_{i=1}^n x_i \\frac1n \\sum_{i=1}^n y_i}{\\frac1n \\sum_{i=1}^n x_i^2 - \\left(\\frac1n \\sum_{i=1}^n x_i\\right)^2}$ \n\net\n\n$b = \\frac1n \\sum_{i=1}^n y_i - \\frac1n \\sum_{i=1}^n x_i a$.\n\n", "metadata": {} }, { "cell_type": "markdown", "source": "Calculons les termes $\\frac1n \\sum_{i=1}^n x_i$, $\\frac1n \\sum_{i=1}^n x_i^2$, $\\frac1n \\sum_{i=1}^n y_i$, $\\frac1n \\sum_{i=1}^n x_i y_i$ :", "metadata": {} }, { "cell_type": "code", "source": "x_mean = np.mean(x_3)\nx2_mean = np.mean(x_3**2)\ny_mean = np.mean(y_train)\nxy_mean = np.mean(x_3 * y_train)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "Vérifions l'unicité de la solution : ", "metadata": {} }, { "cell_type": "code", "source": "x2_mean - x_mean**2", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "Calculons les coefficients de notre droite : ", "metadata": {} }, { "cell_type": "code", "source": "a = (xy_mean - x_mean * y_mean) / (x2_mean - x_mean**2)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "b = y_mean - a * x_mean", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "Nous pouvons maintenant ajouter notre droite de régression au graphique précédent :", "metadata": {} }, { "cell_type": "code", "source": "plt.scatter(x_3, y_train, label=\"Données\")\n\nx_3_min = np.min(x_3) - 5\nx_3_max = np.max(x_3) + 5\nplt.plot([x_3_min, x_3_max], [(a * x_3_min + b), (a * x_3_max + b)], lw=2, c='orange', label=\"Modèle\")\n\nplt.xlabel(\"Longueur de la palette natatoire (mm)\")\nplt.ylabel(\"Poids (g)\")\nplt.legend()", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "### 2 : Régression linéaire multivariée sans scikit-learn (forme close)", "metadata": {} }, { "cell_type": "markdown", "source": "Nous allons maintenant utiliser nos trois variables.\n\nPour simplifier les notations, nous allons ajouter une colonnes de 1 à notre matrice $X$ :\n\n$$\n\\underbrace{\n\\begin{bmatrix}\n x_{11} & x_{12} & \\dots & x_{1p} \\\\\n x_{21} & x_{22} & \\dots & x_{2p} \\\\\n \\vdots & \\vdots & \\ddots & \\vdots \\\\\n x_{n1} & x_{n2} & \\dots & x_{np}\n\\end{bmatrix}}_{X} \\to\n\\underbrace{\n\\begin{bmatrix}\n 1 & x_{11} & x_{12} & \\dots & x_{1p} \\\\\n 1 & x_{21} & x_{22} & \\dots & x_{2p} \\\\\n \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n 1 & x_{n1} & x_{n2} & \\dots & x_{np}\n\\end{bmatrix}}_{\\text{$X$ with bias variable}}\n$$\n\nCe qui nous permet de noter le modèle de la façon suivante :\n$f: \\boldsymbol{x} \\in \\mathbb{R}^{p+1} \\mapsto \\boldsymbol{\\beta}^\\top \\boldsymbol{x}$ with $\\boldsymbol{\\beta} \\in \\mathbb{R}^{p+1}$.", "metadata": {} }, { "cell_type": "code", "source": "X_ones = np.hstack((np.ones((X_train.shape[0], 1)), X_train))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "La minimisation du risque empirique s'écrit alors\n$$\\boldsymbol{\\beta}^* \\in \\arg\\min \\frac{1}{n}\\sum_{i=1}^n \\left(y_i - \\boldsymbol{\\beta}^\\top \\boldsymbol{x}_i \\right)^2 = \\frac{1}{n}\\left(\\boldsymbol{y} - X \\boldsymbol\\beta\\right)^\\top(\\boldsymbol{y} - X\\boldsymbol\\beta)$$\n\nC'est un problème d'optimisation convexe, que l'on résout en annulant le gradient de la fonction à minimiser.\n\nOn obtient alors que la solution $\\boldsymbol\\beta^*$ vérifie\n$(X^\\top X) \\boldsymbol\\beta^* = X^\\top \\boldsymbol{y}.$\n\n__Si $X^\\top X$ est inversible__, on obtient une unique solution\n$\\boldsymbol\\beta^* = (X^\\top X)^{-1} X^\\top \\boldsymbol{y}.$", "metadata": {} }, { "cell_type": "code", "source": "beta_star = np.dot(np.linalg.inv(np.dot(X_ones.T, X_ones)), np.dot(X_ones.T, y_train))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "print(\"Le poids (g) d'un manchot est prédit par %.2f + %.2f x bill_length_mm + %.2f x bill_depth_mm + %.2f x flipper_length_mm\" % (tuple(beta_star)))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "### 3 : Régression logistique sans scikit-learn", "metadata": {} }, { "cell_type": "markdown", "source": "Dans le cas de la régression logistique, la minimisation du risque empirique revient à minimiser\n\n$$J(\\boldsymbol\\beta) = \\frac1n \\sum_{i=1}^n - y_i \\log(\\sigma(\\boldsymbol\\beta^\\top \\boldsymbol{x}_i))- (1-y_i) \\log(1-\\sigma(\\boldsymbol\\beta^\\top \\boldsymbol{x}_i)) $$\n\nComme $\\sigma^\\prime(u) = \\sigma(u) (1-\\sigma(u))$, on obtient comme gradient pour $J$ :\n$$\\nabla J(\\boldsymbol\\beta) = - \\frac1n \\sum_{i=1}^n \\left(y_i - \\sigma(\\boldsymbol\\beta^\\top \\boldsymbol{x}_i) \\right) \\boldsymbol{x}_i$$\n\nCependant, il n'existe pas de solution explcite à \n$$\\frac1n \\sum_{i=1}^n \\left(y_i - \\sigma(\\boldsymbol\\beta^\\top \\boldsymbol{x}_i) \\right) \\boldsymbol{x}_i = 0,$$\nwe do not have a closed-form solution.\n\nNous avons donc recours à la méthode du gradient.\n\nLa __méthode du gradient__ permet d'obtenir une solution en progressant itérativement dans la direction opposée au gradient de la fonction à minimiser. Plus précisément, à chaque itération $t$, le vecteur $\\boldsymbol\\beta$ est mis à jour par :\n$$\\boldsymbol\\beta^{(t+1)} = \\boldsymbol\\beta^{(t)} - \\alpha \\nabla_{\\boldsymbol\\beta} J(\\boldsymbol\\beta^{(t)}).$$\noù $J: \\mathbb{R}^d \\rightarrow \\mathbb{R} $ est la fonction à minimiser.\n\nDans le cas de la régression linéaire, $d=p+1$ et $J(\\boldsymbol\\beta) = \\frac{1}{n}\\left(\\boldsymbol{y} - X \\boldsymbol\\beta\\right)^\\top(\\boldsymbol{y} - X\\boldsymbol\\beta).$\n\n$\\alpha$ est le __pas__ de la descente de gradient ; en apprentissage, on parle de __vitesse d'apprentissage__ (_learning rate_ en anglais).\n\nOn arrête d'itérer :\n* soit quand un nombre prescrit d'itérations est atteint ;\n* soit quand la norme du gradient est moins qu'un seuil prescrit, appelé __tolérance__.", "metadata": {} }, { "cell_type": "markdown", "source": "C'est la méthode que nous allons mettre en œuvre.", "metadata": {} }, { "cell_type": "markdown", "source": "#### Calcul du gradient", "metadata": {} }, { "cell_type": "code", "source": "def sigmoid(u):\n \"\"\"\n Fonction sigmoide.\n \n Paramètres :\n ------------\n u: float\n Nombre réel.\n \"\"\"\n return (1./(1+np.exp(-u)))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "def logistic_loss(X, y, b_vector):\n \"\"\"\n Risque empirique pour une régression logistique.\n \n Paramètres\n ---------- \n X: (n_samples, n_features+1) numpy array\n La matrice de données.\n \n y: (n_samples, ) numpy array \n Le vecteur d'étiquettes\n \n b_vector: (n_features+1, ) numpy array\n Le vecteurs de coefficients de la régression logistique\n \"\"\"\n # Partie du risque correspondant aux étiquettes positives\n where_y_pos = np.where(y==1)[0]\n loss_pos = - np.sum(np.log(sigmoid(X[where_y_pos, :].dot(b_vector))))\n\n # Partie du risque correspondant aux étiquettes négatives \n where_y_neg = np.where(y==0)[0]\n loss_neg = - np.sum(np.log(1 - sigmoid(X[where_y_neg, :].dot(b_vector))))\n \n return (loss_pos + loss_neg)/np.size(y)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "def evaluate_gradient_logistic(X, y, b_vector):\n \"\"\" \n Gradient du risque empirique de la régression logistique\n \n Paramètres\n ---------- \n X: (n_samples, n_features+1) numpy array\n La matrice de données.\n \n y: (n_samples, ) numpy array \n Le vecteur d'étiquettes\n \n b_vector: (n_features+1, ) numpy array\n Le vecteurs de coefficients de la régression logistique\n \"\"\"\n num_samples = X.shape[0]\n diff = sigmoid(X.dot(b_vector)) - y \n return np.sum(np.multiply(diff, X.T), axis=1)/num_samples", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "#### Méthode du gradient", "metadata": {} }, { "cell_type": "code", "source": "def gradient_descent_logistic(X, y, learning_rate=1e-1, max_iters=30, tol=1e-2):\n \"\"\"\n Méthode du gradient pour la régression logistique. \n \n Paramètres\n ---------- \n X: (n_samples, n_features+1) numpy array\n La matrice de données.\n \n y: (n_samples, ) numpy array \n Le vecteur d'étiquettes\n \n learning_rate: float\n La vitesse d'apprentissage\n \n max_iters: int\n Le nombre maximal d'itérations\n \n tol: float\n La tolérance\n \"\"\"\n # Initialisation aléatoire des coefficients de régression\n beta_current = np.random.rand(X.shape[1])\n #beta_current = np.hstack((10*np.random.rand(1), np.random.rand(X.shape[1]-1)))\n \n # Liste pour stocker les valeurs de la fonction de perte\n # à chaque itération\n losses = [logistic_loss(X, y, beta_current)] \n # Liste pour stocker les valeurs des coefficients de régression\n # à chaque itération\n betas = [beta_current.copy()]\n\n for i in range(max_iters): \n # Mise à jour de beta_current\n gradient = evaluate_gradient_logistic(X, y, beta_current) \n beta_current = beta_current - learning_rate * gradient\n \n # Ajouter la fonction de perte à la liste losses\n losses.append(logistic_loss(X, y, beta_current))\n \n # Ajouter la valeur des coefficients de régression à la liste betas\n betas.append(beta_current.copy())\n\n # Vérifier si la tolérance est atteinte\n if np.linalg.norm(gradient) < tol: \n break\n\n return(np.array(losses), betas)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "#### Application aux données", "metadata": {} }, { "cell_type": "code", "source": "(X_train, X_test, y_train, y_test) = model_selection.train_test_split(X, y_classif, test_size=0.3, random_state=25, \n stratify=y_classif)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "X_ones = np.hstack((np.ones((X_train.shape[0], 1)), X_train))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "tolerance = 1e-4\niterations = 35\nlr = 2e-4\nlosses, betas = gradient_descent_logistic(X_ones, y_train, learning_rate=lr, max_iters=iterations, tol=tolerance)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "plt.plot(np.arange(len(losses)), losses, 'o-')\n\nplt.xlabel(\"number of iterations\")\nplt.ylabel(\"value of the loss\")\nplt.title(\"Batch gradient descent\")", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "print(\"La probabilité qu'un manchot soit mâle est prédite par sigma(\", end='')\nprint(\"%.2f + %.2f x bill_length_mm + %.2f x bill_depth_mm + %.2f x flipper_length_mm)\" % (tuple(betas[-1])))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "__Remarque :__ Le risque empirique est infini à l'initialisation, d'où un `RuntimeWarning: divide by zero encountered in log`. En pratique, cette approche pourrait être améliorée par :\n* un pré-traitement des données consistant à centrer-réduire chaque variable ;\n* l'utilisation d'un taux d'apprentissage variable (élevé au début, plus faible à la fin).", "metadata": {} }, { "cell_type": "markdown", "source": "#### Performance", "metadata": {} }, { "cell_type": "code", "source": "X_test_ones = np.hstack((np.ones((X_test.shape[0], 1)), X_test))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "y_test_pred = np.where(sigmoid(X_test_ones.dot(betas[-1])) > 0.5, 1, 0)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "confusion_matrix_display = metrics.ConfusionMatrixDisplay(metrics.confusion_matrix(y_test, y_test_pred))\nconfusion_matrix_display.plot()", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "print(\"%d manchots mâles ont été incorrectement prédits femelle.\" % metrics.confusion_matrix(y_test, y_test_pred)[1, 0])\nprint(\"%d manchots femelles ont été incorrectement prédits mâles.\" % metrics.confusion_matrix(y_test, y_test_pred)[0, 1])", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "print(\"%.f%% des prédictions du modèle sur le jeu de test sont correctes.\" % (100*metrics.accuracy_score(y_test, y_test_pred)))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "__Remarque :__ On confirme ici que la procédure de descente de gradient que nous avons implémentée n'est pas optimale.", "metadata": {} }, { "cell_type": "markdown", "source": "#### Méthode du gradient avec données centrées et réduites", "metadata": {} }, { "cell_type": "markdown", "source": "Nous allons prétraiter les données de manière à ce que chaque variable de notre matrice $\\boldsymbol X$ soit centrée et réduite, c'est à dire que chaque colonne de $\\boldsymbol X$ (à part la première, qui sert de constante) a une moyenne de 0 et un écart-type de 1.\n\nAinsi, l'initialisation de la méthode de gradient fournit des estimées qui sont autour de $0$, c'est-à-dire dans la région où la fonction $\\sigma$ ne sature pas (i.e., son gradient n'est proche de zéro).", "metadata": {} }, { "cell_type": "code", "source": "from sklearn import preprocessing", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "scaler = preprocessing.StandardScaler().fit(X_train)\nX_scaled = scaler.transform(X_train)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "print('means: [%.3f, %.3f, %.3f]' % (tuple(X_scaled.mean(axis = 0))))\nprint('standard deviations: [%.3f, %.3f, %.3f]' % tuple(X_scaled.std(axis = 0)))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "X_ones = np.hstack((np.ones((X_scaled.shape[0], 1)), X_scaled))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "tolerance = 1e-4\niterations = 35\nlr = 2\nlosses, betas = gradient_descent_logistic(X_ones, y_train, learning_rate=lr, max_iters=iterations, tol=tolerance)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "plt.plot(np.arange(len(losses)), losses, 'o-')\n\nplt.xlabel(\"number of iterations\")\nplt.ylabel(\"value of the loss\")\nplt.title(\"Batch gradient descent\")", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "On veut estimer le modèle suivant : $y \\sim \\sigma(\\beta_0 + x_1 \\beta_1 + x_2 \\beta_2 + x_3 \\beta_3)$. Or, puisque nous avons effecuté un changement de variable, nous avons à la place estimé le modèle suivant : $$y \\sim \\sigma(\\tilde{\\beta_0} + u_1 \\tilde{\\beta_1} + u_2 \\tilde{\\beta_2} + u_3 \\tilde{\\beta_3}),$$ où $u_i = \\frac{x_i - \\mu_i}{\\sigma_i}$ sont les nouvelles variables, et $\\mu_i$ et $\\sigma_i$ sont les moyennes et écart-types respectifs de $x_i$.\n\nPar conséquent, il faut effectuer le changement de variable suivant sur les coefficients obtenus :\n$$\\beta_0 = \\tilde{\\beta_0} - \\tilde{\\beta_1}\\frac{\\mu_1}{\\sigma_1} - \\tilde{\\beta_2}\\frac{\\mu_2}{\\sigma_2} - \\tilde{\\beta_3}\\frac{\\mu_3}{\\sigma_3},$$\n$$\\beta_i = \\tilde{\\beta_i} / \\sigma_i, 1\\leq i \\leq 3$$", "metadata": {} }, { "cell_type": "code", "source": "beta_scaled = betas[-1]", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "beta0 = beta_scaled[0] - np.sum(beta_scaled[1:4] * scaler.mean_ / scaler.scale_)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "beta_original = np.insert(beta_scaled[1:4] / scaler.scale_, 0, beta0)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "print(\"La probabilité qu'un manchot soit mâle est prédite par sigma(\", end='')\nprint(\"%.2f + %.2f x bill_length_mm + %.2f x bill_depth_mm + %.2f x flipper_length_mm)\" % (tuple(beta_original)))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "__Remarque__: Les coefficients sont très éloignés de ceux qu'on avait obtenus sans prétraitement des données.", "metadata": {} }, { "cell_type": "markdown", "source": "#### Performance après avoir centré-réduit les données", "metadata": {} }, { "cell_type": "code", "source": "X_test_scaled = np.hstack((np.ones((X_test.shape[0], 1)), X_test))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "y_test_pred = np.where(sigmoid(X_test_scaled.dot(beta_original)) > 0.5, 1, 0)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "confusion_matrix_display = metrics.ConfusionMatrixDisplay(metrics.confusion_matrix(y_test, y_test_pred))\nconfusion_matrix_display.plot()", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "print(\"%d manchots mâles ont été incorrectement prédits femelle.\" % metrics.confusion_matrix(y_test, y_test_pred)[1, 0])\nprint(\"%d manchots femelles ont été incorrectement prédits mâles.\" % metrics.confusion_matrix(y_test, y_test_pred)[0, 1])", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "print(\"%.f%% des prédictions du modèle sur le jeu de test sont correctes.\" % (100*metrics.accuracy_score(y_test, y_test_pred)))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "Nous obtenons un bien meilleur modèle avec un preprocessing des données, qui permet, comme nous l'avons vu, d'éviter les problèmes de convergences que peut rencontrer la méthode de descente de gradient.", "metadata": {} }, { "cell_type": "markdown", "source": "### KNN 1 : Sélection du nombre de plus proches voisins par validation croisée", "metadata": {} }, { "cell_type": "markdown", "source": "Dans ce notebook, nous avons fixé le nombre de plus proches voisins $k$ (`n_neighbors` dans les paramètres de `KNeighborsClassifier` ou `KNeighborsRegressor`) à $k=7$. \n\nNous allons maintenant voir comment __sélectionner__ la meilleure valeur de $k$ dans une liste. On parle ici de __sélection de modèle__.\n\n$k$ est un __hyperparamètre__ du modèle ; l'algorithme des k plus proches voisins nous dit comment procéder à $k$ fixé et le choix de $k$ ne fait pas partie de la procédure. De la même façon, le nombre de couches intermédiaires dans un réseau de neurones est un hyperparamètre. À l'inverse, les coefficients d'une régression linéaire sont des __paramètres__ du modèle, car c'est eux que la procédure d'apprentissage/entraînement/_fit_ cherche à déterminer.\n\nOn peut imaginer approcher cette question en évaluant, sur le jeu de test, la performance de plusieurs $k$NNs, utilisant différentes valeurs de $k$. Il s'agirait alors de choisir la valeur de $k$ donnant la meilleure performance. Cependant, la performance sur le jeu de test du $k$NN sélectionné ne reflète pas bien la performance en généralisation modèle : les données que nous utilisons pour l'évaluation ne sont plus utilisées qu'à ce seul but, puisqu'elles ont aussi été utilisées pour la sélection du modèle.\n\nIl faut donc utiliser uniquement le jeu d'entraînement pour la sélection de modèle.\n\nOn peut de nouveau séparer ce jeu de données en deux jeux, un jeu d'entraînement et un jeu de validation. On évaluera alors les performances de $k$NNs entraînés sur ce nouveau jeu d'entraînement en comparant leurs performances sur le jeu de validation.\n\nCependant, si nous utilisons de nouveau les proportions 70/30 pour créer ce nouveau jeu d'entraînement, on utilise alors 49% des données seulement pour l'entraînement ! C'est peu, surtout si le jeu de données est assez petit ; en effet, plus on a de données, mieux on apprend. On préfère ainsi en pratique utiliser une __validation croisée__ sur le jeu d'entraînement. Il s'agit de découper le jeu d'entraînement en C blocs (appelés _folds_ en anglais ; typiquement C=5 ou C=10), et de faire une série de C expériences dans lesquelles un des blocs est utilisé comme jeu de validation et les (C-1) autres blocs forment le jeu d'entraînement. Ainsi, chaque bloc de données est utilisé une fois pour la validation et (C-1) fois pour l'entraînement. On obtient ainsi C mesures de performance (pour chaque valeur de $k$), que l'on moyenne, et on choisit la valeur de $k$ donnant la meilleure performance moyenne.", "metadata": {} }, { "cell_type": "markdown", "source": "#### GridSearchCV", "metadata": {} }, { "cell_type": "markdown", "source": "Cette procédure est automatisée dans [la classe `GridSearchCV`](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html) du module `model_selection` de `scikit-learn`", "metadata": {} }, { "cell_type": "markdown", "source": "On procède alors de la façon suivante :", "metadata": {} }, { "cell_type": "markdown", "source": "1. Fixons une liste de valeurs à évaluer pour l'hyperparamètre $k$\n\nNous choisissons ici uniquement des valeurs impaires, pour faciliter la détermination de la classe majoritaire dans le cas d'une classification binaire (avec un nombre pair de voisin, on peut ne pas avoir de classe majoritaire).", "metadata": {} }, { "cell_type": "code", "source": "n_neighbors_values = [5, 7, 9, 13, 15, 17, 19, 21]", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "2. Initialisons un objet de la classe `GridSearchCV`, auquel nous passons comme paramètres :\n* la classe de modèle qui nous intéresse, ici `KNeighborsClassifier` pour le problème de classification ;\n* un dictionnaire indiquant, pour chacun des paramètres de ce modèle à évaluer (ici nous en avons un seul, `n_neighbors`), la liste de valeurs à tester (`n_neighbors_values`) ;\n* le nombre de _folds_ (ici 5) ;\n* la fonction de score à utiliser pour déterminer la performance d'un modèle ; ici nous allons utiliser l'_accuracy_. Pour plus de possibilités, voir [la documentation](https://scikit-learn.org/stable/modules/model_evaluation.html#the-scoring-parameter-defining-model-evaluation-rules).", "metadata": {} }, { "cell_type": "code", "source": "knnclass_cv = model_selection.GridSearchCV(neighbors.KNeighborsClassifier(),\n {'n_neighbors': n_neighbors_values}, cv=5, scoring='accuracy')", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "3. Faisons tourner la validation croisée sur le jeu d'entraînement", "metadata": {} }, { "cell_type": "code", "source": "knnclass_cv.fit(X_train, y_train)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "#### Nombre optimal de voisins ", "metadata": {} }, { "cell_type": "markdown", "source": "Nous avons maintenant accès au nombre optimal de voisins ainsi qu'à l'_accuracy_ correspondante, qui sont des attributs de `knnclass_cv` : ", "metadata": {} }, { "cell_type": "code", "source": "print(\"La valeur optimale de k est %d, pour une accuracy cross-validée de %.f%%.\" % (knnclass_cv.best_params_['n_neighbors'], (100*knnclass_cv.best_score_)))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "#### Performance du modèle optimal", "metadata": {} }, { "cell_type": "markdown", "source": "Nous avons aussi accès au modèle optimal, c'est-à-dire le modèle entraîné sur l'intégralité du jeu d'entraînement, avec la valeur optimale de $k$, dans `knnclass_cv.best_estimator_`. Nous pouvons donc évaluer sa performance sur le jeu de test.", "metadata": {} }, { "cell_type": "code", "source": "y_test_pred = knnclass_cv.best_estimator_.predict(X_test)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "metrics.plot_confusion_matrix(knnclass_cv.best_estimator_, X_test, y_test)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "print(\"%d manchots mâles ont été incorrectement prédits femelle.\" % metrics.confusion_matrix(y_test, y_test_pred)[1, 0])\nprint(\"%d manchots femelles ont été incorrectement prédits mâles.\" % metrics.confusion_matrix(y_test, y_test_pred)[0, 1])", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "print(\"%.f%% des prédictions du modèle sur le jeu de test sont correctes.\" % (100*metrics.accuracy_score(y_test, y_test_pred)))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "### KNN 2 : Plus proche voisin (au singulier) sans scikit-learn", "metadata": {} }, { "cell_type": "markdown", "source": "Nous allons maintenant voir comment implémenter l'algorithme du plus proche voisin ($k$=1) sans utiliser `scikit-learn`.", "metadata": {} }, { "cell_type": "code", "source": "def my_nearest_neighbor(x, Xtr, ytr):\n \"\"\"\n Retourne l'étiquette du plus proche voisin dans Xtr de x.\n \n Parameters\n ----------\n x: np.array of shape (p, )\n L'observation à étiquetter\n \n Xtr: np.array of shape (n, p)\n La matrice de données d'entrainement\n \n ytr: np.array of shape (n, )\n Les étiquettes des données d'entrainement. \n \"\"\"\n \n # Calculer les carrés distances euclidiennes entre x et chaque ligne de Xtr\n squared_distances = [np.sum((x - Xtr[idx, :])**2) for idx in range(Xtr.shape[0])]\n \n # Trouver l'indice de la plus petite de ces valeurs \n nn_idx = np.argmin(squared_distances)\n \n # Retourner l'étiquette correspondante\n return ytr[nn_idx]", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "Nous pouvons maintenant utiliser cette fonction sur toutes les observations de `X_test` pour obtenir les prédictions :", "metadata": {} }, { "cell_type": "code", "source": "my_y_test_pred = [my_nearest_neighbor(X_test[idx, :], X_train, y_train) for idx in range(X_test.shape[0])]", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "Nous pouvons comparer ces prédictions à celles de `scikit-learn` :", "metadata": {} }, { "cell_type": "code", "source": "onnclass = neighbors.KNeighborsClassifier(n_neighbors=1)\nonnclass.fit(X_train, y_train)\ny_test_pred = onnclass.predict(X_test)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "my_y_test_pred - y_test_pred", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "### KNN 3 : Plus proches voisins (au pluriel) sans scikit-learn", "metadata": {} }, { "cell_type": "markdown", "source": "#### Pour la classification", "metadata": {} }, { "cell_type": "code", "source": "(X_train, X_test, y_train, y_test) = model_selection.train_test_split(X, y_classif, test_size=0.3, random_state=25)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "def my_k_nearest_neighbors_classif(x, Xtr, ytr, k):\n \"\"\"\n Retourne l'étiquette de la plus fréquente parmi les k plus proches voisins dans Xtr de x.\n \n Parameters\n ----------\n x: np.array of shape (p, )\n L'observation à étiquetter\n \n Xtr: np.array of shape (n, p)\n La matrice de données d'entrainement\n \n ytr: np.array of shape (n, )\n Les étiquettes des données d'entrainement. \n \n k: int\n Nombre de plus proches voisins\n \"\"\"\n \n # Calculer les carrés distances euclidiennes entre x et chaque ligne de Xtr\n squared_distances = [np.sum((x - Xtr[idx, :])**2) for idx in range(Xtr.shape[0])]\n \n # Trouver les indices des k plus petites valeurs dans cette liste,\n # en récupérant les k premiers indices de la liste triée\n closest_indices = np.argsort(squared_distances)[:k]\n \n # Retourner l'étiquette la plus fréquente dans cette liste\n # qui est 1 si la moyenne des étiquettes est supérieure à 0.5 et 0 sinon\n return (1 if (np.mean(ytr[closest_indices])) > 0.5 else 0)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "Nous pouvons maintenant utiliser cette fonction sur toutes les observations de `X_test` pour obtenir les prédictions :", "metadata": {} }, { "cell_type": "code", "source": "my_y_test_pred = [my_k_nearest_neighbors_classif(X_test[idx, :], X_train, y_train, 7) for idx in range(X_test.shape[0])]", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "Nous pouvons comparer ces prédictions à celles de `scikit-learn` :", "metadata": {} }, { "cell_type": "code", "source": "knnclass = neighbors.KNeighborsClassifier(n_neighbors=7)\nknnclass.fit(X_train, y_train)\ny_test_pred = knnclass.predict(X_test)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "y_test_pred = knnclass.predict(X_test)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "my_y_test_pred - y_test_pred", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "#### Pour la régression", "metadata": {} }, { "cell_type": "code", "source": "(X_train, X_test, y_train, y_test) = model_selection.train_test_split(X, y_regress, test_size=0.3, random_state=25)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "def my_k_nearest_neighbors_regress(x, Xtr, ytr, k):\n \"\"\"\n Retourne la moyenne des étiquettes des k plus proches voisins dans Xtr de x.\n \n Parameters\n ----------\n x: np.array of shape (p, )\n L'observation à étiquetter\n \n Xtr: np.array of shape (n, p)\n La matrice de données d'entrainement\n \n ytr: np.array of shape (n, )\n Les étiquettes des données d'entrainement. \n \n k: int\n Nombre de plus proches voisins\n \"\"\"\n \n # Calculer les carrés distances euclidiennes entre x et chaque ligne de Xtr\n squared_distances = [np.sum((x - Xtr[idx, :])**2) for idx in range(Xtr.shape[0])]\n \n # Trouver les indices des k plus petites valeurs dans cette liste,\n # en récupérant les k premiers indices de la liste triée\n closest_indices = np.argsort(squared_distances)[:k]\n \n #print(list(closest_indices))\n \n # Retourner la moyenne des étiquettes dans cette liste\n return (np.mean(ytr[tuple([closest_indices])]))", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "Nous pouvons maintenant utiliser cette fonction sur toutes les observations de `X_test` pour obtenir les prédictions :\n\nAttention, `y_train` est une série pandas ; pour obtenir l'array numpy correspondant, nous utilisons `y_train.values`.", "metadata": {} }, { "cell_type": "code", "source": "my_y_test_pred = [my_k_nearest_neighbors_regress(X_test[idx, :], X_train, y_train.values, 7) for idx in range(X_test.shape[0])]", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": "Nous pouvons comparer ces prédictions à celles de `scikit-learn` :", "metadata": {} }, { "cell_type": "code", "source": "knnreg = neighbors.KNeighborsRegressor(n_neighbors=7)\nknnreg.fit(X_train, y_train)\ny_test_pred = knnreg.predict(X_test)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "y_test_pred = knnreg.predict(X_test)", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "code", "source": "my_y_test_pred - y_test_pred", "metadata": {}, "execution_count": null, "outputs": [] } ] }