{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*M2 MIGS -- Statistique pour les big data -- Université de Bourgogne -- 2022-2023*\n",
    "\n",
    "# Travaux pratiques 2 : Chemin des solutions de l'estimateur LASSO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Ces travaux pratiques sont notés. Durée : 3h, tous les documents autorisés, usage d'internet autorisé (sauf pour communiquer avec un étudiant ou une autre perssonne).** A la fin de la séance il faudra envoyer ce document avec pour titre : \"TP_Big_Data_NOM_Prénom\" à Patrick.Tardivel@u-bourgogne.fr  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# on importe les modules numpy et pyplot\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from math import sqrt\n",
    "import PIL\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Dans le cellule ci-dessous ajouter votre code permettant de résoudre le problème LASSO :**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercice 1 \n",
    "\n",
    "On sonsidère la matrice $X=\\begin{pmatrix} 4 & 2 & 1\\end{pmatrix}$, $y=4$ et $\\lambda=2$. \n",
    "\n",
    "**1** Déterminer le polyèdre nul du LASSO $\\{r\\in \\mathbb{R}: \\|X^Tr\\|_\\infty\\}$ et donner le résidu $\\hat r$ du LASSO.\n",
    "\n",
    "**2** En déduire la solution du LASSO $$\\min_{b\\in \\mathbb{R}^3}\\left\\{ \\frac{1}{2}\\|y-Xb\\|_2^2+\\lambda\\|b\\|_1\\right\\}$$ à la main."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Répondre aux questions 1 et 2 dans cette cellule."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3** Vérifier la solution obtenue avec votre méthode de résolution numérique numérique"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution numérique (si votre ne code ne fonctionne contactez-moi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercice 2 : Comparing Ridge and LASSO estimators for data analysis\n",
    "\n",
    "### Première partie : Analyse de l'article\n",
    "\n",
    "L'article  \"Comparing Ridge and LASSO estimators for data analysis\" est un acte de conférence ; c'est-à-dire que ce court papier est un support écrit qui accompagne une présentation à la conférence internationale :\n",
    "\"Information Technology and Nanotechnology\". Quelques questions liées à cet article sont listées ci-dessous :\n",
    "\n",
    "**1)** On considère le modèle de régression linéaire ${\\bf Y}={\\bf XB}+\\mathcal{E}$. Rappeler l'expression explicite de l'estimateur des moindres carrés. Quelle hypothèse est faîte sur ${\\bf X}$ pour que cet estimateur soit bien défini ? D'après les auteurs quel défaut peut avoir l'estimateur des moindres carrés ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2)** Soit $\\lambda> 0$, l'estimateur ridge est défini par \n",
    "$$\\widetilde{\\bf B}_\\lambda :=\\underset{B\\in \\mathbb{R}^k}{\\rm argmin}\\,\\|{\\bf V}-{\\bf W}B\\|_2^2+\\lambda \\|B\\|_2^2.$$\n",
    "L'expression explicite de l'estimateur ridge est donnée dans ce papier. Rappeler cette expression.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3)** D'après les auteurs, quelle est la différence principale entre l'estimateur ridge et l'estimateur LASSO ? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Deuxième partie : chemin des solution du LASSO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le but des questions suivantes est de construire le chemin des solutions du LASSO (il me semble que le chemin des solutions du LASSO donné à la figure 1 de l'article est faux). On rappelle que  le chemin des solutions du LASSO est la fonction \n",
    "$\\lambda>0\\mapsto \\widetilde{\\bf B}_\\lambda$ où $\\widetilde{\\bf B}_\\lambda$ est l'unique solution du problème LASSO :\n",
    "$$\\widetilde{\\bf B}_\\lambda :=\\underset{B\\in \\mathbb{R}^k}{\\rm argmin}\\,\\frac{1}{2}\\|{\\bf V}-{\\bf W}B\\|_2^2+\\lambda \\|B\\|_1.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>fixed acidity</th>\n",
       "      <th>volatile acidity</th>\n",
       "      <th>citric acid</th>\n",
       "      <th>residual sugar</th>\n",
       "      <th>chlorides</th>\n",
       "      <th>free sulfur dioxide</th>\n",
       "      <th>total sulfur dioxide</th>\n",
       "      <th>density</th>\n",
       "      <th>pH</th>\n",
       "      <th>sulphates</th>\n",
       "      <th>alcohol</th>\n",
       "      <th>quality</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>7.4</td>\n",
       "      <td>0.700</td>\n",
       "      <td>0.00</td>\n",
       "      <td>1.9</td>\n",
       "      <td>0.076</td>\n",
       "      <td>11.0</td>\n",
       "      <td>34.0</td>\n",
       "      <td>0.99780</td>\n",
       "      <td>3.51</td>\n",
       "      <td>0.56</td>\n",
       "      <td>9.4</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>7.8</td>\n",
       "      <td>0.880</td>\n",
       "      <td>0.00</td>\n",
       "      <td>2.6</td>\n",
       "      <td>0.098</td>\n",
       "      <td>25.0</td>\n",
       "      <td>67.0</td>\n",
       "      <td>0.99680</td>\n",
       "      <td>3.20</td>\n",
       "      <td>0.68</td>\n",
       "      <td>9.8</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>7.8</td>\n",
       "      <td>0.760</td>\n",
       "      <td>0.04</td>\n",
       "      <td>2.3</td>\n",
       "      <td>0.092</td>\n",
       "      <td>15.0</td>\n",
       "      <td>54.0</td>\n",
       "      <td>0.99700</td>\n",
       "      <td>3.26</td>\n",
       "      <td>0.65</td>\n",
       "      <td>9.8</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>11.2</td>\n",
       "      <td>0.280</td>\n",
       "      <td>0.56</td>\n",
       "      <td>1.9</td>\n",
       "      <td>0.075</td>\n",
       "      <td>17.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>0.99800</td>\n",
       "      <td>3.16</td>\n",
       "      <td>0.58</td>\n",
       "      <td>9.8</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>7.4</td>\n",
       "      <td>0.700</td>\n",
       "      <td>0.00</td>\n",
       "      <td>1.9</td>\n",
       "      <td>0.076</td>\n",
       "      <td>11.0</td>\n",
       "      <td>34.0</td>\n",
       "      <td>0.99780</td>\n",
       "      <td>3.51</td>\n",
       "      <td>0.56</td>\n",
       "      <td>9.4</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1594</th>\n",
       "      <td>6.2</td>\n",
       "      <td>0.600</td>\n",
       "      <td>0.08</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.090</td>\n",
       "      <td>32.0</td>\n",
       "      <td>44.0</td>\n",
       "      <td>0.99490</td>\n",
       "      <td>3.45</td>\n",
       "      <td>0.58</td>\n",
       "      <td>10.5</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1595</th>\n",
       "      <td>5.9</td>\n",
       "      <td>0.550</td>\n",
       "      <td>0.10</td>\n",
       "      <td>2.2</td>\n",
       "      <td>0.062</td>\n",
       "      <td>39.0</td>\n",
       "      <td>51.0</td>\n",
       "      <td>0.99512</td>\n",
       "      <td>3.52</td>\n",
       "      <td>0.76</td>\n",
       "      <td>11.2</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1596</th>\n",
       "      <td>6.3</td>\n",
       "      <td>0.510</td>\n",
       "      <td>0.13</td>\n",
       "      <td>2.3</td>\n",
       "      <td>0.076</td>\n",
       "      <td>29.0</td>\n",
       "      <td>40.0</td>\n",
       "      <td>0.99574</td>\n",
       "      <td>3.42</td>\n",
       "      <td>0.75</td>\n",
       "      <td>11.0</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1597</th>\n",
       "      <td>5.9</td>\n",
       "      <td>0.645</td>\n",
       "      <td>0.12</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.075</td>\n",
       "      <td>32.0</td>\n",
       "      <td>44.0</td>\n",
       "      <td>0.99547</td>\n",
       "      <td>3.57</td>\n",
       "      <td>0.71</td>\n",
       "      <td>10.2</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1598</th>\n",
       "      <td>6.0</td>\n",
       "      <td>0.310</td>\n",
       "      <td>0.47</td>\n",
       "      <td>3.6</td>\n",
       "      <td>0.067</td>\n",
       "      <td>18.0</td>\n",
       "      <td>42.0</td>\n",
       "      <td>0.99549</td>\n",
       "      <td>3.39</td>\n",
       "      <td>0.66</td>\n",
       "      <td>11.0</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1599 rows × 12 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      fixed acidity  volatile acidity  citric acid  residual sugar  chlorides  \\\n",
       "0               7.4             0.700         0.00             1.9      0.076   \n",
       "1               7.8             0.880         0.00             2.6      0.098   \n",
       "2               7.8             0.760         0.04             2.3      0.092   \n",
       "3              11.2             0.280         0.56             1.9      0.075   \n",
       "4               7.4             0.700         0.00             1.9      0.076   \n",
       "...             ...               ...          ...             ...        ...   \n",
       "1594            6.2             0.600         0.08             2.0      0.090   \n",
       "1595            5.9             0.550         0.10             2.2      0.062   \n",
       "1596            6.3             0.510         0.13             2.3      0.076   \n",
       "1597            5.9             0.645         0.12             2.0      0.075   \n",
       "1598            6.0             0.310         0.47             3.6      0.067   \n",
       "\n",
       "      free sulfur dioxide  total sulfur dioxide  density    pH  sulphates  \\\n",
       "0                    11.0                  34.0  0.99780  3.51       0.56   \n",
       "1                    25.0                  67.0  0.99680  3.20       0.68   \n",
       "2                    15.0                  54.0  0.99700  3.26       0.65   \n",
       "3                    17.0                  60.0  0.99800  3.16       0.58   \n",
       "4                    11.0                  34.0  0.99780  3.51       0.56   \n",
       "...                   ...                   ...      ...   ...        ...   \n",
       "1594                 32.0                  44.0  0.99490  3.45       0.58   \n",
       "1595                 39.0                  51.0  0.99512  3.52       0.76   \n",
       "1596                 29.0                  40.0  0.99574  3.42       0.75   \n",
       "1597                 32.0                  44.0  0.99547  3.57       0.71   \n",
       "1598                 18.0                  42.0  0.99549  3.39       0.66   \n",
       "\n",
       "      alcohol  quality  \n",
       "0         9.4        5  \n",
       "1         9.8        5  \n",
       "2         9.8        5  \n",
       "3         9.8        6  \n",
       "4         9.4        5  \n",
       "...       ...      ...  \n",
       "1594     10.5        5  \n",
       "1595     11.2        6  \n",
       "1596     11.0        6  \n",
       "1597     10.2        5  \n",
       "1598     11.0        6  \n",
       "\n",
       "[1599 rows x 12 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "wine = pd.read_csv(\"winequality-red.csv\", sep=';')\n",
    "wine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standardisation des données\n",
    "\n",
    "y = wine.to_numpy()[:,11]\n",
    "y_mean = np.mean(y)\n",
    "y_std = np.std(y)\n",
    "V = (y-y_mean)/(sqrt(1599)*y_std)\n",
    "\n",
    "X = wine.to_numpy()[:,0:11]\n",
    "X_mean = np.mean(X, axis=0)\n",
    "X_std = np.std(X, axis=0)\n",
    "W = np.hstack((np.ones((1599,1)),(X-X_mean)/(sqrt(1599)*X_std)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**4)** Calculer $\\lambda_{\\max}=\\|{\\bf W}^T{\\bf V}\\|_{\\infty}$. Que peut-on dire de l'estimateur LASSO lorsque $\\lambda \\ge \\lambda_{\\max}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5)** On pose $\\widetilde B_\\lambda$ l'unique solution du problème LASSO pour le paramètre de régularisation $\\lambda$. Construire le chemin des solutions de l'estimateur LASSO $\\widetilde B_\\lambda$, c'est-à-dire tracer sur un même graphique les courbes des douze composantes de cet estimateur :\n",
    "    $$ \\lambda>0\\mapsto [\\widetilde B_\\lambda]_i \\text{ pour }i=0, \\dots, 11.$$\n",
    "Pour ce graphique on discrétisera l'intervalle $[ 0,001\\lambda_\\max ; 1,2\\lambda_\\max]$ en 1000 points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**6)** Sur le même graphique  ajouter les points de coordonnées $(0,\\widehat {\\bf B}_i)$ où $\\widehat {\\bf B}$ est l'estimateur des moindres carrés de ${\\bf B}$ calculé avec ${\\bf W}$ et ${\\bf V}$. Que remarquez-vous ? Justifier votre commentaire. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Troisième partie : choix du paramètre de régularisation pour minimiser l'erreur de prédiction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'erreur de prédiction de l'estimateur LASSO est $\\mathbb{E}(\\|W\\widetilde B_\\lambda - \\mathbb{E}(V)\\|^2_2)$. La formule SURE (Stein Unbiased Risk Estimate), facilement calculable, fournit un estimateur sans biais de l'erreur de prédiction du LASSO : $\\|V-W\\widetilde B_\\lambda\\|_2^2 - n \\sigma^2 + 2\\sigma^2 \\|\\widetilde B_\\lambda\\|_0$. Afin d'appliquer cette formule on doit dans un premier temps estimer $\\sigma^2$. \n",
    "\n",
    "**7)** A partir de l'estimateur des moindres carrés, donner une estimation de $\\sigma^2$. \n",
    "\n",
    "**8)** Quel paramètre de régularisation $\\lambda$ appartenant discrétisation de l'intervalle $[ 0,001\\lambda_\\max ; 1,2\\lambda_\\max]$\n",
    " minimise la formule SURE ? Pour ce paramètre de régularisation, quelles sont les variables explicatives non sélectionnées par l'estimateur LASSO ?  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
