{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TP 6 : La Régression Linéaire"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le but du TP est de mettre en oeuvre une régression linéaire sur un exemple jouet puis sur des données réelles. Pour cela, vous aurez besoin des formules du cours pour retrouver les coefficients permettant de calculer une régression linéaire. \n",
    "\n",
    "\n",
    "### Compétences associées :\n",
    "\n",
    "| Numéro        | Compétence           | \n",
    "|:------------- |:--------------------:|\n",
    "|RL001 | Savoir définir un modèle linéaire simple|\n",
    "|RL004 | Être capable de modéliser la prédiction d’une nouvelle donnée|\n",
    "|RL005 | Être capable de calculer l’erreur de prédiction|\t\n",
    "|RL007|Savoir calculer un intervalle de prédiction de la régression linéaire|\n",
    "|RL103|Savoir retrouver l’expression des coefficients optimaux en fonction de X et y|\n",
    "|RL104|Savoir retrouver l’expression de la prédiction d’une donnée x en fonction des coefficients optimaux|\n",
    "|RL201|Savoir calculer le R² en fonction de la décomposition de la variance\n",
    "|RL202|Savoir interpréter le R²|\n",
    "|PY402|Être capable de calculer les coefficients optimaux de la régression linéaire (avec numpy)|\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex. 1 : La régression simple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Nous donnons les couples d’observations suivantes :\n",
    "\n",
    "x| 18| 7 | 14| 31| 21|  5| 11| 16| 26| 29\n",
    "-|---|---|---|---|---|---|---|---|---|---\n",
    "y| 55| 17| 36| 85| 62| 18| 33| 41| 63| 87\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x=np.array([18, 7, 14, 31, 21, 5, 11, 16, 26, 29, 8])\n",
    "y=np.array([55, 17, 36, 85, 62, 18, 33, 41, 63, 87, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Tracer le graphique des couples $x$ et $y$. A partir de ce graphe, peut on soupconner une relation linéaire entre les variables $x$ et $y$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "plt.plot(x,y,'o')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Déterminer pour ces observations la droite de régression au sens des moindres carrés. Pour cela il faudra déterminer les coefficients a et b."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "sx2 = np.sum((x-np.mean(x))**2)\n",
    "sxy=((x-np.mean(x)).T @ (y-np.mean(y)))\n",
    "a = sxy/sx2\n",
    "b = np.mean(y) - a*np.mean(x)\n",
    "n = len(x)\n",
    "print(a,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Tracez la droite de régression sur les données."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "xm = np.min(x)-6; \n",
    "xM = np.max(x)+6;\n",
    "ym = a*xm+b;\n",
    "yM = a*xM+b;\n",
    "plt.plot([xm,xM],[ym,yM],'g-', linewidth=2) # trace une ligne entre [xm, ym] et[xM,yM]\n",
    "plt.plot(x,y,'o')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Donner une estimation des erreurs $\\epsilon_i,i=1,n$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "yp = a*x+b\n",
    "e = y - yp\n",
    "e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "6. Calculer la moyenne empirique et la variance empirique des erreurs. Comparez la variance avec l'estimateur non biaisé de la variance donné par $s_y^2 = \\frac{1}{n-1}\\sum_{i=1}^n (y_i - \\overline y)^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "n = len(e)\n",
    "sigma2_hat = 1/(n-1)* e.T@e # la moyenne des résidus est égale à 0\n",
    "print(sigma2_hat) #estimateur non biaisé de la variance (si ils l'ont vu)\n",
    "print(np.var(e))\n",
    "print(np.mean(e))\n",
    "plt.hist(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "7. Donner une estimation plausible de $y$ lorsque $x = 17$. Quel intervalle de confiance associer à cette prédiction ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "from scipy.stats import t\n",
    "\n",
    "xp = 17\n",
    "yp = a*xp+b\n",
    "\n",
    "#Les aider la dessus\n",
    "t_s = t.ppf(0.025,n-2)  # statistique de student à n-2 degres de liberté  \n",
    "interval = t_s*np.sqrt(sigma2_hat*(1+1/n+(xp-np.mean(x))**2/sx2));\n",
    "plt.plot([xp, xp],[yp-interval, yp+interval],'c')\n",
    "xm = np.min(x)\n",
    "xM = np.max(x)\n",
    "ym = a*xm+b\n",
    "yM = a*xM+b\n",
    "\n",
    "plt.plot([xm,xM],[ym,yM],'gv-', linewidth=2)\n",
    "plt.plot(x,y,'b+')\n",
    "plt.show()\n",
    "print(interval)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "8. a) Donner une estimation plausible de y lorsque x = 48. Quel intervalle de confiance associer à cette prédiction ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "xp = 48\n",
    "yp = a*xp+b;\n",
    "t_s = t.ppf(0.025,n-2)  # loi de student à n-2 ddl\n",
    "interval = t_s*np.sqrt(sigma2_hat*(1+1/n+(xp-np.mean(x))**2/sx2));\n",
    "print(interval)\n",
    "print(yp+interval, yp-interval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "xp =  np.arange(xm,xM,.1)\n",
    "yp = a*xp+b\n",
    "tv = t.ppf(0.025,n-2)  \n",
    "print(tv)\n",
    "interval = tv*np.sqrt(sigma2_hat*(1+1/n+(xp-np.mean(x))**2/sx2));\n",
    "plt.plot([xp, xp],[yp-interval, yp+interval],'c')\n",
    "plt.plot([xm,xM],[ym,yM],'gv-', linewidth=2)\n",
    "plt.plot(x,y,'b+')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "plt.plot(np.abs(interval))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b) Comparez le avec le précédent\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "source": [
    "On s'écarte de la moyenne de x, donc ça augmente"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "9. Une nouvelle observation nous est fournie : $x_{11} = 48$ et $y_{11} = 2$\n",
    "    1. Que devient la droite de regression linéaire ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "x_new = np.append(x, 48)\n",
    "y_new= np.append(y, 2)\n",
    "sx2 = np.sum((x_new-np.mean(x_new))**2),\n",
    "sy2 = np.sum((y_new-np.mean(y_new))**2),\n",
    "sxy=((x_new-np.mean(x_new)).T @ (y_new-np.mean(y_new)))\n",
    "a = sxy/sx2\n",
    "b = np.mean(y) - a*np.mean(x)\n",
    "xm = np.min(x_new)-1;\n",
    "xM = np.max(x_new)+1;\n",
    "ym = a*xm+b;\n",
    "yM = a*xM+b;\n",
    "plt.plot(x_new,y_new,'b+')\n",
    "plt.plot([xm, xM],[ym, yM],'r');\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "B. Quelle est l’influence de ce point sur cette droite ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "source": [
    "Il influe de manière disproportionné sur la droite de régression.\n",
    "Calcul des $w_i$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "w = (x_new - np.mean(x_new)) / np.sum((x_new - np.mean(x_new))**2)\n",
    "print(w)\n",
    "plt.bar(np.arange(len(w)), w);\n",
    "#on voit bien la forte influence du dernier point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "## Ex 2 : Régression sur $CO_2$ et température moyenne du globe\n",
    "\n",
    "Le but de ce second exercice est d'appliquer la régression au sens des moindres carrés sur des données réelles. Ici la variable explicative est le taux de $CO_2$ dans l'atmosphère, et la variable à expliquer la température moyenne à la surface du globe, avec pour référence  à 0 la température moyenne sur la période 1961-1990. \n",
    "Les données sont les mêmes que dans le TP4, où nous avons trouvé une corrélation d'environ 0.8. \n",
    "Maintenant, nous allons chercher à calculer les estimateurs du modèle linéaire. Il faut donc appliquer tout ce que nous avons vu sur des données jouets, et interprétez ces résultats.\n",
    "\n",
    "1. Chargez les données contenues dans 'data_global_warming.npz'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "data = np.load('data_global_warming.npz')\n",
    "temp = data['temp']\n",
    "co2 = data['co2']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "plt.plot(co2,temp,'o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Créez une fonction `regression_fit` calculant les estimateurs $a$ et $b$ donné $x$ et $y$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def regression_fit(x_,y_):\n",
    "    ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "#Essayez de leur faire comprendre qu'il faut être efficace en termes de complexité\n",
    "def regression_fit(x_,y_):\n",
    "    import numpy as np\n",
    "    mx = np.mean(x_)\n",
    "    my = np.mean(y_)\n",
    "    sx2 = np.sum((x_- mx)**2)\n",
    "    sxy= np.sum((x_- mx) * (y_-my))\n",
    "    a = sxy/sx2\n",
    "    b = my - a*mx\n",
    "    return a,b\n",
    "a,b = regression_fit(co2,temp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Créez une fonction `regression_predict` qui calcule l'estimation de $y$ donnés $x$, $a$ et $b$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def regression_predict(x_,a_,b_):\n",
    "    ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "def regression_predict(x_,a_,b_):\n",
    "    return x_*a_ + b_\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Créez une fonction `residuals` qui calcule les résidus. Quelles sont les paramètres de cette fonction ? Est ce que les erreurs vous paraissent avoir une distribution gaussienne ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def residuals(y_true, y_pred):\n",
    "    ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "def residuals(y_true, y_pred):\n",
    "    return y_true-y_pred\n",
    "\n",
    "erreurs = residuals(temp,regression_predict(co2,a,b))\n",
    "plt.hist(erreurs)\n",
    "print(np.mean(erreurs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Calculez les valeurs prédites sur l'intervalle \\[1,5\\] avec un pas de 0.1 et plottez le résultat sur les données."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "x_plot = np.arange(1,5.1,.1)\n",
    "y_plot = regression_predict(x_plot,a,b)\n",
    "plt.plot(co2,temp,'o')\n",
    "plt.plot(x_plot,y_plot)\n",
    "print(y_plot.shape,x_plot.shape,b.shape,b.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Affichez la valeur de $a$ et concluez sur la nature de la relation entre $CO_2$ et température."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "source": [
    "a est > 0 donc la temp augmente avec le taux de CO2. \n",
    "\n",
    "En conclusion, essayez de les questionner si un modele lineaire est le plus adapté ? Surement qu'un modele polynomial serait mieux. Les interroger sur comment faire avec un modèle linéaire. La réponse : projeter x vers $[x,x^2, x^3, ..., ]$ ou même avec des combinaisons, et avppliquer une régression linéaire \n",
    "multiple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bonus: Différentes manières de calculer le $R^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = regression_predict(co2,a,b)\n",
    "\n",
    "(np.sum((y_pred - np.mean(temp))**2) /np.sum((temp - np.mean(temp))**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1 - (np.sum((temp - y_pred)**2)/np.sum((temp-np.mean(temp))**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.corrcoef(y_pred.T, temp.T)[0,1]**2"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "iml",
   "language": "python",
   "name": "iml"
  },
  "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}