{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TP 7 - Régression Linéaire Multiple\n",
    "\n",
    "\n",
    "### Compétences associées :\n",
    "\n",
    "| Numéro        | Compétence           | \n",
    "|:------------- |:--------------------:|\n",
    "|RL101|\tSavoir définir un modelé linéaire multiple|\n",
    "|RL102| Savoir écrire la représentation matricielle du problème de régression au sens des moindres carrés|\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",
    "|PY403| Être capable de calculer les coefficients optimaux de la régression linéaire  (avec scikit-learn)|\t\t\n",
    "\n",
    "## Ex 1: Régression multiple\n",
    "\n",
    "Nous allons évaluer notre modèle de régression linéaire sur le jeu de données diabetes. Nous allons utiliser la variable à expliquer correspondant à un indice d'évolution de la maladie.\n",
    "Plus de détails directement sur le site de sklearn : [lien](https://scikit-learn.org/stable/datasets/toy_dataset.html#diabetes-dataset)\n",
    "\n",
    "Pour implémenter la régression multiple, nous allons utiliser la modélisation matricielle.\n",
    "\n",
    "1. Charger les données et faites les analyses statistiques et visualisations qui vous paraissent appropriées. Commentez la deuxième colonne des données."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "from sklearn.datasets import load_diabetes\n",
    "..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "from sklearn.datasets import load_diabetes\n",
    "X,y = load_diabetes(return_X_y=True)\n",
    "n,p = X.shape\n",
    "plt.boxplot(X); # la deuxieme colonne est une variable qualitative (sexe) qui a été quantifiée (deux valeurs)\n",
    "plt.figure();\n",
    "plt.hist(y);\n",
    "\n",
    "plt.figure()\n",
    "cov = X.T@X\n",
    "plt.matshow(cov)\n",
    "plt.colorbar()\n",
    "print(np.mean(X,axis=0))\n",
    "print(np.var(X,axis=0)*n) # -> cf la doc pour l'explication"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Construisez la matrice $X$ pour la régression linéaire. Pour rappel, dans le formalisme matriciel, il faut ajouter une colonne de 1 afin de prendre en compte le coefficient $a_0$.\n",
    "\n",
    "La fusion de tableau sur numpy peut paraitre compliquée, étant donné le nombre de façons de faire différentes. La liste des fonctions existantes est disponible [ici](https://numpy.org/doc/stable/reference/routines.array-manipulation.html#joining-arrays). Un tutoriel est dispo [ici](https://kanoki.org/2020/01/10/concatenating-arrays-in-numpy/).\n",
    "Pour cet exercice, on utilisera la fonction `np.concatenate` pour fusionner `X` et une vecteur colonne de 1 obtenu par `np.ones((n,1))`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "n,p = X.shape\n",
    "X_reg = np.concatenate((X,np.ones((n,1))), axis=1)\n",
    "#équivalent à \n",
    "X_reg2 = np.column_stack((X,np.ones((n,1))))\n",
    "print(X_reg.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Calculez le vecteur $\\alpha$ optimisant le problème de régression au sens des moindres carrés. Visualisez les différentes valeurs avec `plt.bar`\n",
    " \n",
    "Attention : Il est généralement déconseillé d'inverser une matrice si cela n'est pas absolument nécessaire. Vous devriez trouver la fonction qu'il vous faut [ici](https://numpy.org/doc/stable/reference/routines.linalg.html#solving-equations-and-inverting-matrices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "alpha = np.linalg.solve(X_reg.T@X_reg, X_reg.T@y) #Attention à bien utiliser X_reg !!!\n",
    "plt.bar(np.arange(p+1),alpha)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Pour chaque observation, calculez la prédiction faite par le modèle en calculant $z=X\\alpha$. \n",
    "Comparez aux valeurs réelles, et calculez la moyenne des valeurs absolues des résidus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "corrrection"
    ]
   },
   "outputs": [],
   "source": [
    "z = X_reg@alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "erreurs = z-y\n",
    "np.mean(z-y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. Calculez le coefficient de détermination R² pour évaluer la qualité de votre prédiction. Commentez votre résultat. Quelles pistes pour l'améliorer ?\n",
    "\n",
    "Bonus : Vous pouvez essayer de le calculer de plusieurs manières différentes"
   ]
  },
  {
   "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((z-np.mean(y))**2)\n",
    "print(SCE/SCT)\n",
    "print(np.corrcoef(z,y)[0,1]**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "# comparer avec sklearn\n",
    "from sklearn.linear_model import LinearRegression\n",
    "linreg = LinearRegression()\n",
    "linreg.fit(X,y)\n",
    "z=linreg.predict(X)\n",
    "print(linreg.coef_,linreg.intercept_)\n",
    "print(alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex. 2 — Histoire de température\n",
    "Un méteorologue nous a fourni les données suivantes de température et de concentration en ozone mesurées au centre de la ville de Caen, à la même heure, les lundis de printemps 2007."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|x : Température | y : Ozone|\n",
    "|:---:| --- |\n",
    "|23.1 | 81.4|\n",
    "|18.9 | 52.6|\n",
    "|16.4 | 48.4|\n",
    "|14.6 | 38.0|\n",
    "|20.3 | 39.4|\n",
    "|25.4 | 34.2|\n",
    "|21.5 | 61.9|\n",
    "|22.1 | 77.3|\n",
    "|17.3 | 34.0|\n",
    "|14.8 | 49.8|\n",
    "|25.5 | 97.0|\n",
    "|23.1 | 67.8|\n",
    "|24.2 | 92.7|"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "pour expliquer la concentration d’ozone en fonction de la température, les météorologues proposent le modèle suivant :\n",
    "\\begin{equation}\n",
    "y = \\begin{cases}\n",
    "    b + \\varepsilon & \\text{ si } x \\lt c\\\\\n",
    "    ax + d + \\varepsilon & \\text{ si } x \\ge  c\n",
    "  \\end{cases}\n",
    "\\end{equation}\n",
    "\n",
    "Ces données sont présentes dans le fichier meteo.csv sur Moodle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Réécrire le modèle en éliminant d en utilisant la contrainte de continuité du modèle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "data = np.loadtxt(\"meteo.csv\", delimiter=\",\")\n",
    "x = data[:, 0]\n",
    "print(x)\n",
    "y = data[:, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "source": [
    "La continuité implique $b = ac+d$ soit $d = b - ac$. Sans $d$ le modèle s'écrit alors : \n",
    "\\begin{equation}\n",
    "y = \\begin{cases}\n",
    "    b + \\varepsilon & \\text{ si } x \\lt c\\\\\n",
    "    a(x-c) + b + \\varepsilon & \\text{ si } x \\ge  c\n",
    "  \\end{cases}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Sachant que c = 21, Nous allons calculer les coefficients a et b optimaux au sens des moindres carrés.\n",
    "\n",
    "    a) Expliciter  la transformation de X tel que le modèle s’écrive $y = X \\alpha + \\varepsilon$. Transformez vos données pour appliquer cette transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "c=21\n",
    "x_transform = data[:, 0].copy()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "source": [
    "x_transform = 0 si x < c, x-c sinon.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "# on aurait pu faire en compréhension de liste, mais ce n'est pas le but de M8\n",
    "for i,xi in enumerate(x_transform):\n",
    "    if (xi < c) :\n",
    "        x_transform[i] = 0\n",
    "    else :\n",
    "        x_transform[i] = xi - c\n",
    "\n",
    "plt.plot(x,y,'o')\n",
    "plt.figure()\n",
    "plt.plot(x_transform,y,'o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b) Calculer l’estimateur de a et b au sens des moindres carrés"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "\n",
    "X = np.stack([x_transform, np.ones(len(x_transform))], axis=-1)\n",
    "a = np.linalg.solve(X.T@X, X.T@y)\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "c) Visualiser les observations et le modèle de régression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt \n",
    "plt.plot(x,y,'or')\n",
    "z = X@a\n",
    "plt.plot(x,z,'og')\n",
    "\n",
    "xp = np.array([14,21,28]) # les points charnières\n",
    "# pour les deux premiers, le modeles juste f(x) = b\n",
    "yp = np.array([a[1],a[1],a[0]*(xp[2]-c)+a[1]])\n",
    "\n",
    "plt.plot(xp,yp,'b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3\\. Prédictions\n",
    "    \n",
    " a) Quelle serait d’après votre modèle la concentration d’ozone lorque la température est de 27 degrés."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_x = 27"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "yp = np.array([new_x-c,1]) @ a\n",
    "yp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b) Quelle serait d’après votre modèle la concentration d’ozone lorque la température est de 13 degrés."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "# 13 < 21 -> pas de x dans le modele\n",
    "yp =  np.array([0,1]) @ a\n",
    "yp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4\\. Quelle est la qualité de la régression ? Calculez le coefficient de détermination $R^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "e = y-X@a\n",
    "SCT = np.sum((y-np.mean(y))**2)\n",
    "SCM = np.sum((X*a-np.mean(y))**2)\n",
    "SCR = e.T @ e;\n",
    "R2 = 1 - SCR/SCT\n",
    "SCT, SCM, SCR, R2\n",
    "print(R2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus : Ex. 3 — Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cherchant à expliquer la température d’un four (y) en fonction de l’hygrométrie (x) un boulanger a réalisé dix huit expériences différentes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "x = np.array([-0.20, 0.72, 0.08, 1.39, 0.22, 0.87, 1.38, 0.94, 0.36, 1.22, 0.74, 0.13, 0.21, 0.56, -0.35, -0.30, 0.48,\n",
    "0.48])\n",
    "y = np.array([-0.04, 1.30, 0.08, 1.37, 0.27, 1.36, 1.44, 1.28, 0.62, 1.29, 1.41, 0.12, 0.26, 1.16, -0.04, -0.16, 0.93,\n",
    "0.83])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Le voisin du boulanger, lui même boulanger à la retraite, explique que le modèle linéaire n’est pas le bon. Le bon modèle est le suivant :\n",
    "\n",
    "\\begin{equation}\n",
    "y = \\begin{cases}\n",
    "    ax^2 + bx + \\epsilon  & \\text{ si } x \\le 0.6\\\\\n",
    "    c + \\epsilon & \\text{ si } x \\gt  0.6\n",
    "  \\end{cases}\n",
    "\\end{equation}\n",
    "    a) Ecrire le modèle sous la forme\n",
    "\n",
    "\\begin{equation}\n",
    "        y = X\\alpha + \\epsilon\n",
    "\\end{equation}\n",
    "et construire la matrice X associée"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "source": [
    "La contrainte de continuité implique $a 0.36 + b 0.6 = c$. En éliminant $c$  le modèle devient : \n",
    "\\begin{equation}\n",
    "y = \\begin{cases}\n",
    "    ax^2 + bx + \\epsilon  & \\text{ si } x \\le 0.6\\\\\n",
    "    0.36 a + 0.6 b + \\epsilon & \\text{ si } x \\gt  0.6\n",
    "  \\end{cases}\n",
    "\\end{equation}\n",
    "et donc\n",
    "$$\n",
    "X = \\left(\n",
    " \\begin{array}{cc}\n",
    "  .2^2 & -.2 \\\\\n",
    " .36 & .6 \\\\\n",
    " .08^2 & .08 \\\\\n",
    " .36 & .6  \\\\\n",
    "\\vdots & \\vdots\\\\\n",
    " .48^2 & .48 \n",
    " \\end{array}\n",
    " \\right) \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "ip = np.where(x>0.6)[0]\n",
    "X = np.stack([x**2,x], axis=1)\n",
    "X[ip] = np.ones((len(ip),1))*[0.36, 0.6]\n",
    "print(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2\\. Estimer les paramètres du modèle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "a = np.linalg.solve(X.T@X, X.T@y)\n",
    "e = y - X@a\n",
    "S2hat = e.T@e/16"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3\\. Visualiser les observations et le modèle 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,y,'or')\n",
    "z = X@a\n",
    "plt.plot(x,z,'og')\n",
    "\n",
    "t = np.arange(-0.4, 0.61, 0.01)\n",
    "tp = np.arange(0.6, 1.4, 0.01)\n",
    "\n",
    "t1 = np.stack([t**2, t])\n",
    "tp1 = np.stack([np.ones((len(tp)))*0.36, np.ones((len(tp)))*0.6])\n",
    "\n",
    "Xt = np.concatenate([t1, tp1], axis=1)\n",
    "\n",
    "plt.plot(np.concatenate([t, tp]), Xt.T @ a,'b')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}