{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sélection de variables pour la régression\n",
    "\n",
    "Le but de ce TP est d'implémenter divers algorithmes pour sélectionner uniquement les variables les plus importantes pour un problème de régression. Ces méthodes servent pour améliorer la performance d'un modèle tout en réduisant sa complexité\n",
    "\n",
    "| Numéro        | Compétence           | \n",
    "|:------------- |:--------------------:|\n",
    "|RL206|\tComprendre la définition du Cp de Mallows|\n",
    "|RL207|\tEtre capable de définir un algorithme de sélection de modèle donné un critère |\n",
    "|PY503|\tSavoir implémenter un algorithme itératif de sélection de variables|\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ex. 1 — De l’importance des variables\n",
    "\n",
    "Récupérez le fichier cereal.mat sur moodle. Ce fichier contient les mesures des ingrédients de céréales pour petit déjeuner. Sur chaque type de céréale, 8 variables ont été mesurées. Ces variables sont (dans l’ordre) :\n",
    "- Calories: la teneur en calories (le nombre pour un repas)\n",
    "- Carbo: la teneur en carbohydrates (en gramme pour un repas)\n",
    "- Cups: le nombre de tasse recommandé par repas\n",
    "- Fat: la teneur en graisse (en gramme pour un repas)\n",
    "- Fiber: la teneur en fibres alimentaires (en gramme pour un repas)\n",
    "- Potass: la teneur en potassium (en mg pour un repas)\n",
    "- Protein: la teneur en proteine (en gramme pour un repas)\n",
    "- Sugars: la teneur en sucre (en gramme pour un repas)\n",
    "    \n",
    "Ce fichier peut se décomposer en une matrice `X` qui contient toutes les données, et le tableau `Name` qui contient les noms des marques.\n",
    "\n",
    "À chaque type de céréale a été associé son nom contenu dans la variable Name. Ainsi la première céréale est mat[\"name\"][0] -> ’100% Bran’."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1) Afin de déterminer quelle est la variable la plus importante pour la régression linéaire, effectuer 7 régressions linéaire simple de chacune des autres variables contre  La variable FAT.\n",
    "\n",
    "On commence par charger les données du fichier \"cereal.mat\" et construire la matrice `X` du modèle complet et le vecteur `y` de la variable à expliquer. Vous pouvez jeter un oeil au fichier à l'entete du fichier `cereal.mat`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.io as sio\n",
    "import numpy as np\n",
    "mat = sio.loadmat(\"cereal.mat\")\n",
    "X = mat[\"X\"]\n",
    "y = X[:, 3] # la quatrieme colonne\n",
    "X = np.delete(X, 3, 1)\n",
    "np.unique(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b. Constuire une matrice A permettrant de stocker les paramètres de tous les modèles, ainsi qu'un vecteur $s2$ et un vecteur $R2$ permettant de stocker les $S^2$ et $R^2$ de tous les modèles que l'on souhaite tester.\n",
    "\n",
    "**aide** : `np.zeros`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "A = np.zeros((p,2))\n",
    "s2 = np.zeros((p,))\n",
    "R2 = np.zeros((p,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = ...\n",
    "s2 = ...\n",
    "R2 = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "c. En utilisant la fonction `ma_reg` du précédent TP, calculez les paramètres, les $S^2$ (variance non biaisée) et les $R^2$ (coefficient de détermination) de tous les modèles. \n",
    "\n",
    "Retrouvez quelle est la variable la plus importante au sens du $R^2$\n",
    "\n",
    "**aide** `np.argmax` renvoie l'indice du max plutot que la valeur du max"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(p):\n",
    "    Xn = ...\n",
    "    A[i,:], s2[i], R2[i], diagd = ...\n",
    "\n",
    "with np.printoptions(precision=3, suppress=True):\n",
    "    print('Erreur résiduelle :',s2) \n",
    "    print('               R2 :',R2)\n",
    "\n",
    "v = ...\n",
    "ind = ...\n",
    "\n",
    "print(f\"Le plus grand R2 est {v:3.3f}\")\n",
    "print(f\"La variable la plus importante est la variable {ind+1}\")\n"
   ]
  },
  {
   "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",
    "    import numpy as np\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-p)\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "for i in range(p):\n",
    "    Xn = np.stack([X[:,i], np.ones((n,))], axis=-1)\n",
    "    # calculer la matrice des variables explicatives\n",
    "    A[i,:], s2[i], R2[i], diagd = ma_reg(Xn,y)#  faire la regression\n",
    "\n",
    "with np.printoptions(precision=3, suppress=True):\n",
    "    print('Erreur résiduelle :',s2) \n",
    "    print('               R2 :',R2)\n",
    "\n",
    "v = np.max(R2)\n",
    "ind = np.argmax(R2)\n",
    "\n",
    "print(f\"Le plus grand R2 est {v:3.3f}\")\n",
    "print(f\"La variable la plus importante est la variable {ind+1}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Nous allons maintenant calculer l’estimation du risque statistique (le Cp de Mallows) associé à ce modèle.\n",
    "    \n",
    "    a. calculez la variance estimée sur le modèle complet  (Attention aux dégrès de libertés)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "X_reg= np.hstack([X, np.ones((n,1))])\n",
    "a_complet,s2_complet, R2_complet, _ = ma_reg(X_reg,y)\n",
    "\n",
    "z = X_reg @ a_complet\n",
    "e = y-z\n",
    "#e.T@e  = \\sum e_i ² \n",
    "s2_r = ((e.T@e)/(n-p-1))\n",
    "print(s2_complet,s2_r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s2_r = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b. calculez le $C_p$ du meilleur modèle à une seule variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "p0 = 2\n",
    "X_best_1 = np.hstack([X[:,ind].reshape(-1,1), np.ones((n,1))])\n",
    "A_best_1 = A[ind,:]\n",
    "z_0 = X_best_1 @ A_best_1\n",
    "e_0 = y-z_0\n",
    "Cp_1variable = (1/s2_r) * (e_0.T@e_0) - n + 2 * p0\n",
    "print(Cp_1variable)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Cp_1variable = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Nous allons maintenant trouver quelles sont les deux variables optimales au sens du $C_p$ \n",
    "\n",
    "    a. combien de modèles possible avec deux variables ?  Attention, l'ordre des variables n'a pas d'importance. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "nb_modeles = (p*(p-1))/2\n",
    "nb_modeles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_modeles = ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b. Calculez le $C_p$ de chacun de ces modèles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "s2_2_var = np.ones((p,p))*np.inf\n",
    "r2_2_var = np.ones((p,p))*np.inf\n",
    "for i in range(p):\n",
    "    for j in range(i+1,p):\n",
    "        Xn = np.stack([X[:,i],X[:,j], np.ones((n,))], axis=-1)\n",
    "        _, s2_2_var[i,j], r2_2_var[i,j], diagd = ma_reg(Xn,y)\n",
    "\n",
    "with np.printoptions(precision=2,suppress=True):\n",
    "    print(s2_2_var)\n",
    "    print(r2_2_var)\n",
    "    \n",
    "p0 = 3    \n",
    "Cp_2variable = ((s2_2_var*(n-p0))/s2_r) - n + 2 * p0\n",
    "with np.printoptions(precision=2,suppress=True):\n",
    "    print(Cp_2variable)\n",
    "    \n",
    "\n",
    "best_ind_2 = np.argmin(Cp_2variable)\n",
    "best_cp_2 = Cp_2variable[np.unravel_index(best_ind_2, Cp_2variable.shape)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Estimer le meilleur modèle au sens du $C_p$ par *forward* sélection. \n",
    "\n",
    "N'hésitez pas à réfléchir à l'algorithme sur une feuille de papier auparavant. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "# vide -> on rajoute la meilleure variable \n",
    "p=7\n",
    "var_in = [] # variables selectionnées\n",
    "var_out = [x for x in range(p)] # variables candidates \n",
    "R2_all_models = []\n",
    "Cp = [] # critere final pour le meilleur modele\n",
    "\n",
    "#Calcul du s2 modele complet\n",
    "X_complet = np.c_[X, np.ones(y.shape)]\n",
    "_, s2_complet, _, _ = ma_reg(X_complet,y) \n",
    "\n",
    "# p tours de boucle pour ajouter p variables au fur et à mesure\n",
    "for i in range(p):\n",
    "    # parcours des variables restantes à ajouter\n",
    "    s2_i = []\n",
    "    R2_i = []\n",
    "    for j in range(len(var_out)): # test des variables candidates\n",
    "        #Creation du jeu de données avec X[:,var_out[j]] en plus de X[:,var_in] \n",
    "        liste_variables_courante = var_in + [var_out[j]]\n",
    "        X_courant = X[:,liste_variables_courante]\n",
    "        #Ajout des 1\n",
    "        X_courant = np.concatenate([X_courant, np.ones((len(X_courant), 1))], axis=1)\n",
    "        \n",
    "        #Calcul de la regression, sauvegarde des résultats\n",
    "        _, s2_courant, R2_courant, _ = ma_reg(X_courant,y) \n",
    "        s2_i.append(s2_courant)\n",
    "        R2_i.append(R2_courant)\n",
    "\n",
    "    #On récupère le meilleur modele (s2, R2 et Cp sont équivalents pour un même nombre de variables)\n",
    "    v = np.min(s2_i)\n",
    "    ind  = np.argmin(s2_i)\n",
    "    R2_all_models.append(R2_i[ind])\n",
    "    pi = X_courant.shape[1] # \n",
    "    #Calcul et sauvegarde du Cp pour choix final\n",
    "    Cp.append((n-pi)*v/s2_complet - n + 2*pi)\n",
    "    print(f\"Variable gagnante pour le tour {i} : {var_out[ind]} avec un R2 de {R2_i[ind]} et un Cp de {Cp[-1]}\")\n",
    "    # Mise à jour de var_in et var_out\n",
    "    var_in = var_in + [var_out[ind]]\n",
    "    var_out.remove(var_out[ind])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "correction"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(range(p),Cp)\n",
    "#plt.xlabel = var_in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5\\. **Bonus** Estimer le meilleur modèle au sens du $C_p$ par *backward* sélection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "venvm8-ioPrAr7d-py3.10",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}