{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# TP 8 : Diagnostic de la régression : analyse des observations\n", "Le but du TP est de réaliser un diagnostic d'un modèle de régression d'un point de vue général, ainsi que via l'observation des variables. Ce TP vous permettra d'implémenter le calcul du $R^2$, ainsi que des coefficients de levier et de contribution.\n", "\n", "| Numéro | Compétence | \n", "|:------------- |:--------------------:|\n", "|PY502|\tAfficher un résumé du diagnostic de la régression et savoir l’analyser|\n", "|RL203|\tSavoir analyser les résidus de régression|\n", "|RL204|\tComprendre la définition de l’effet levier d’une observation|\n", "|RL205|\tComprendre la définition de la contribution d’un couple|\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ex. 1 — Détection de points aberrants dans la régression\n", "On considère les n = 11 observations suivantes :\n", "\n", "i| 1| 2 |3| 4| 5| 6| 7| 8| 9 |10| 11\n", "---|---|---|---|---|---|---|---|---|---|---|---\n", "x | -1.0782 |0.0838| 0.1524 |0.2290| 0.4427| 0.4505| 0.5383| 0.8258| 0.9133| 0.9961| 3.0597\n", "y| -0.8307| 0.1538 |0.0926| 0.2336| 0.3903| 0.1005| 0.4812 |0.6595| 0.7175| 0.7648 |0.3519" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les observations $y$ sont liées aux $x$ à travers le modèle linéaire :\n", "$y_i = ax_i + b + \\epsilon_i \\text{, avec } i = 1, n$\n", "où les $\\epsilon_i$ sont des variables aléatoires i.i.d. distribuées selon une loi normale centrée de variance\n", "$\\sigma^2$ inconnue et $(a, b)$ sont deux paramètres inconnus.\n", "\n", "1. Effectuons les calculs liés au modèle de régression et à son diagnostic.\n", "\n", " a) construire une matrice $X$ et un vecteur $y$ contenant les données de la régression que l’on souhaite effectuer." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "data = np.array([\n", " [1, -1.0782, -.8307],\n", " [2, 0.0838, 0.1538],\n", " [3, 0.1524, 0.0926],\n", " [4, 0.2290, 0.2336],\n", " [5, 0.4427, 0.3903],\n", " [6, 0.4505, 0.1005],\n", " [7, 0.5383, 0.4812],\n", " [8, 0.8258, 0.6595],\n", " [9, 0.9133, 0.7175],\n", " [10, 0.9961, 0.7648],\n", " [11, 3.0597, 0.3519]\n", "])\n", "X = ...\n", "y = ...\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "n, _ = data.shape\n", "p = 1\n", "X=np.stack([np.ones((n,)), data[:,1]], axis=-1)\n", "y = data[:,2]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "a) Calculer l’estimation de $a$ et $b$ sens des moindres carrés" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "...\n", "a = ...\n", "b = ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "alpha = np.linalg.solve(X.T@X,X.T@y)\n", "a = alpha[1]\n", "b = alpha[0]\n", "a, b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "b) Plotter $x$ et $y$, ainsi que votre droite de régression." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "plt.plot(X[:,1],y,'o')\n", "z = X@alpha\n", "plt.plot(x,z,'x-')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "c) Calculer des résidus et une estimation non biaisée de la variance des erreurs $\\sigma^2$.\n", "> **Aide** Combien de degré de libertés pour les résidus ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "e = y - z\n", "s2 = np.sum(e**2)/(n-2) #Moyenne des residus à 0; n -p-1 ddl, ici p=1\n", "print(s2)\n", "print(np.mean(e))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "d) Évaluer la qualité globale de la régression" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "SCT = np.sum((y-np.mean(y))**2)\n", "SCE = np.sum(e**2)\n", "R2 = 1 - SCE/SCT #plein de manieres de calculer\n", "R2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "e) Calculer le levier de chaque observation. Affichez les valeurs sous forme d'un graphique en barres." ] }, { "cell_type": "markdown", "metadata": { "tags": [ "correction" ] }, "source": [ "D'après le cours : \n", "$$\n", "H = X (X^\\top X)^{-1}X^\\top\n", "$$\n", "les leviers de chaque point sont les $H_{ii}$\n", " on peut aussi écrire : \n", " $$\n", " H = X G\n", " $$\n", " avec \n", " $$\n", " G = (X^\\top X)^{-1}X^\\top\n", " $$\n", " La matrice $G$ est donc aussi solution des $n$ systèmes linéaires : \n", " $$\n", " (X^\\top X) \\; G = X^\\top\n", " $$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "h = ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "H = X@np.linalg.inv(X.T@X)@X.T\n", "H = X@np.linalg.solve(X.T@X,X.T) # plus stable et plus rapide\n", "h = np.diag(H)\n", "print(h)\n", "plt.bar(range(n),h)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "g) Calculer la contribution de chaque observation. Idem, affichez les valeur sous forme de graphiques en barres." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "c = ..." ] }, { "cell_type": "markdown", "metadata": { "tags": [ "correction" ] }, "source": [ "d'après le cours \n", "$$\n", " c_i = \\frac{H_{ii}}{(1 - H_{ii})^2}\\frac{\\hat \\varepsilon_i^2}{p s^2}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "c = h/((1-h)**2)*(e**2)/(p*s2)\n", "plt.bar(np.arange(n),c)\n", "print(c)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "h) Visualisez tous ces résultats dans un tableau contenant le vecteur des observations\n", "$x$, le vecteur des réponses $y$, les prédiction du modèle, les résidus, les leviers, les\n", "contributions et une variable binaire repérant les points dont la contribution dépasse\n", "$4/n$. Ces résultats vous paraissent ils cohérents ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "res = np.array([x, y, z, e, h, c , c>(4/n)])\n", "print('-------------------------------------------------------')\n", "print(' x y z e h c Aberrant ?')\n", "print('-------------------------------------------------------')\n", "with np.printoptions(precision=3, suppress=True):\n", " print(res.T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. Nous allons maintenant « consolider » ces résultats en écrivant une fonction réutilisable (et par la même nous abstraire du copier/coller et de ces effets néfastes).\n", " \n", "a) Écrire une fonction Python ma_reg permettant d’intégrer les calculs fait précédemment. Votre fonction prendra la matrice $X$ et le vecteur $y$ comme paramètres d’entrée et retournera :\n", " \n", "- l’estimation des paramètres de la régression sens des moindres carrés,\n", "- une estimation de la variance des erreurs $\\sigma^2$,\n", "- le coefficient de détermination $R^2$\n", "- une matrice n × 3 contenants des éléments permettant de réaliser le diagnostic de la régression à savoir :\n", " - les résidus,\n", " - le levier de chaque observation,\n", " - la contribution de chaque observation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def ma_reg(X,y):\n", " pass" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "def ma_reg(X,y): #a voir avec les eleves si ils veulent avoir X avec la colonne de 1 intégré ou pas\n", " '''\n", " X : matrice de taille n par p+1, avec derniere colonne de 1\n", " '''\n", " n,p = X.shape\n", " a = np.linalg.solve(X.T@X,X.T@y)\n", " z = X@a\n", " e = y - z \n", " s2 = np.sum(e**2)/(n-p) # ici le nb de col de X est p+1\n", " SCT = np.sum((y-np.mean(y))**2)\n", " SCE = np.sum(e**2)\n", " R2 = 1 - SCE/SCT\n", " H = X@np.linalg.solve(X.T@X,X.T) # plus stable et plus rapide\n", " h = np.diag(H)\n", " c = h/((1-h)**2)*(e**2)/(p*s2)\n", " \n", " dv = np.stack([e, h, c], axis=-1)\n", " \n", " return a, s2, R2, dv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "c) Écrire un script Python permettant de tester cette fonction. Que constatez vous ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "a, s2, R2, diagd = ma_reg(X,y)\n", "a,s2,R2, diagd\n", "# le dernier point est suspect !" ] }, { "cell_type": "markdown", "metadata": { "tags": [ "correction" ] }, "source": [ "Sans le dernier point, la regression est bien meilleure" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "d) Que se passe t’il si l’on exécute les instruction suivantes ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "yp = X\n", "a, s2, R2, diagd = ma_reg(X,yp)" ] }, { "cell_type": "markdown", "metadata": { "tags": [ "correction" ] }, "source": [ "la fonction ma_reg attend une matrice X de taille $n \\times p$ ce qui est bon, et un vecteur de taille $n$. Ici, à la place de ce vecteur, on passe en paramètre une matrice yp de taille $n \\times p$ puisqu'il s'agit de $X$ !" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "e) Ajouter au début de la fonction ma_reg un test permettant de vérifier que les dimensions de la matrice $X$ et du vecteur $y$ sont bien conformes" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def ma_reg(X,y):\n", " pass" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "def ma_reg(X,y):\n", " '''\n", " X : numpy array de dimention n x p\n", " y : numpy array de dimention n\n", " '''\n", " n,p = X.shape\n", " m = np.prod(y.shape)\n", " if n != m:\n", " raise Exception('X doit etre une matrice de n lignes et p colonnes et y un vecteur de n lignes')\n", " \n", " a = np.linalg.solve(X.T@X,X.T@y)\n", " z = X@a\n", " e = y - z # z - y\n", " s2 = np.sum(e**2)/(n-2)\n", " SCT = np.sum((y-np.mean(y))**2)\n", " SCE = np.sum(e**2)\n", " R2 = 1 - SCE/SCT\n", " H = X@np.linalg.solve(X.T@X,X.T) # plus stable et plus rapide\n", " h = np.diag(H)\n", " c = h/((1-h)**2)*(e**2)/(p*s2)\n", " \n", " dv = np.stack([e, h, c])\n", " \n", " return a, s2, R2, dv\n", "\n", "a, s2, R2, diagd = ma_reg(X,y)\n", "#yp = X\n", "#a, s2, R2, diagd = ma_reg(X,yp)\n", "print(diagd.T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3\\. Effectuer le diagnostic de la régression\n", " \n", "a) trouver une observation avec un résidu important et un faible levier" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { "tags": [ "correction" ] }, "source": [ "C'est le cas de l'avant dernier point. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "b) trouver une observation avec un faible résidu et un fort levier\n", " \n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { "tags": [ "correction" ] }, "source": [ " point 1 \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "c) trouver une observation avec un fort résidu et un levier important. Vérifier que vous avez bien détecté qu’il s’agit d’un point aberrant.\n", " \n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { "tags": [ "correction" ] }, "source": [ " dernier point" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "d) Après avoir éliminé le point aberrant, refaire toutes les étapes de la régression" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "Xn = np.delete(X, -1, 0)\n", "yn = np.delete(y, -1 ,0)\n", "a2, s22, R22, D2 = ma_reg(Xn,yn)\n", "\n", "plt.subplot(2,1,1)\n", "plt.plot(Xn[:,1],yn,'o')\n", "plt.plot(Xn[:,1],Xn@a,'r--')\n", "plt.plot(Xn[:,1],Xn@a2,'g-')\n", "\n", "plt.subplot(2,1,2)\n", "plt.plot(Xn[:,1], D2[0,:],'*m')\n", "plt.grid()\n", "plt.show()\n", "\n", "xa = Xn@a2\n", "print('-------------------------------------------- ')\n", "print(' x y z e h c ')\n", "print('-------------------------------------------- ')\n", "R = np.stack([Xn[:,1], yn[:], xa[:], D2[0,:], D2[1,:], D2[2,:]],axis=-1)\n", "with np.printoptions(precision=3, suppress=True):\n", " print(R)\n", " \n", "plt.plot(x,y,'o')\n", "z = X@a\n", "plt.plot(x,z,'x-')\n", "zn = X@an\n", "plt.plot(x,zn,'x-')\n", "plt.legend(['les observations','la régression acvec tous ','la régression sans le dernier point'])\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "a) trouver une observation avec un résidu important et un faible levier\n", " \n", "C'est le cas du 6ème point dans la régression \n", " \n", "b) trouver une observation avec un faible résidu et un fort levier\n", " \n", "C'est le cas du premier point ($H_{11} = 0.75$)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ex. 2 — Modèle mal spécifié\n", "On considère les n = 11 observations suivantes :\n", "\n", "x| 1.00 |1.50| 2.00| 2.50| 3.00| 3.50| 4.00| 4.50| 5.00| 5.50| 6.0000\n", "---|---|---|---|---|---|---|---|---|---|---|---|\n", "y| 5.10| 4.56| 4.85| 7.48| 9.10| 11.65| 14.87| 19.20| 24.23| 32.10| 39.80\n", "\n", "1. Effectuer le régression linéaire entre x et y et déterminer le coefficient de corrélation et visualiser les résidus. Ces résidus sont-ils conformes aux hypothèses de la régression ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "data = np.array([[1, 5.1],\n", " [1.5, 4.56],\n", " [2, 4.85],\n", " [2.5, 7.48],\n", " [3, 9.1],\n", " [3.5, 11.65],\n", " [4, 14.87],\n", " [4.5, 19.2],\n", " [5, 24.23],\n", " [5.5, 32.1],\n", " [6, 39.8]])\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "\n", "x = data[:,0]\n", "y = data[:,1]\n", "n = len(x)\n", "X = np.stack([np.ones((n,)), x], axis=-1)\n", "\n", "(a, s2, R2i, diagd) = ma_reg(X,y)\n", "z = X@a\n", "e = diagd[0,:]\n", "\n", "plt.subplot(3,1,1)\n", "plt.plot(x,y,'o')\n", "plt.plot(x,z,'r-')\n", "plt.subplot(3,1,2)\n", "plt.plot(x,e,'o')\n", "plt.show() # non, ils ne suivent pas une loi normale mais on une structure" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2\\. Effectuer la régression polynomiale (degré 2) entre x et y et déterminer le coefficient de détermination. La régression polynomiale consiste à faire l’hypothèse que les observations $y$ sont liées aux $x$ à travers le modèle quadratique :\n", "\n", "$y_i = a_0 + a_1x_i + a_2x^2_i + \\epsilon_i\\text{, avec } i = 1, n$\n", "\n", "où les $\\epsilon_i$ sont des variables aléatoires i.i.d. distribuées selon une loi normale centrée de\n", "variance $\\sigma^2$ inconnue et $(a_0, a_1, a_2)$ sont trois paramètres inconnus." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "jupyter": { "outputs_hidden": true }, "tags": [ "correction" ] }, "outputs": [], "source": [ "X = np.stack([np.ones((n,)), x, x**2], axis=-1)\n", "a, s2, R2f, diagd = ma_reg(X,y)\n", "print(R2f)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3\\. Visualiser les résidus de la régression et vérifier qu’ils vérifient bien les hypothèses de la régression." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "plt.subplot(3,1,1)\n", "plt.plot(x,y,'o')\n", "z = X@a\n", "plt.plot(x,z,'r-')\n", "plt.subplot(3,1,2)\n", "e = diagd[0,:]\n", "plt.plot(x,e,'o')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "4\\. Comparez les coefficients de détermination des deux régressions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "correction" ] }, "outputs": [], "source": [ "print(R2i)\n", "print(R2f)" ] } ], "metadata": { "celltoolbar": "Tags", "kernelspec": { "display_name": "mlenv", "language": "python", "name": "mlenv" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.6" } }, "nbformat": 4, "nbformat_minor": 4 }