.. _tdnote2017rst: ================================ 1A.e - TD noté, 16 décembre 2016 ================================ .. only:: html **Links:** :download:`notebook `, :downloadlink:`html `, :download:`python `, :downloadlink:`slides `, :githublink:`GitHub|_doc/notebooks/exams/td_note_2017.ipynb|*` Régression linéaire avec des variables catégorielles. .. code:: ipython3 from jyquickhelper import add_notebook_menu add_notebook_menu() .. contents:: :local: Exercice 1 ---------- On suppose qu’on dispose d’un ensemble d’observations :math:`(X_i, Y_i)` avec :math:`X_i, Y_i \in \mathbb{R}`. La régression linéaire consiste une relation linéaire :math:`Y_i = a X_i + b + \epsilon_i` qui minimise la variance du bruit. On pose : .. math:: E(a, b) = \sum_i (Y_i - (a X_i + b))^2 On cherche :math:`a, b` tels que : .. math:: a^*, b^* = \arg \min E(a, b) = \arg \min \sum_i (Y_i - (a X_i + b))^2 La fonction est dérivable et on trouve : .. math:: \frac{\partial E(a,b)}{\partial a} = - 2 \sum_i X_i ( Y_i - (a X_i + b)) \text{ et } \frac{\partial E(a,b)}{\partial b} = - 2 \sum_i ( Y_i - (a X_i + b)) Il suffit alors d’annuler les dérivées. On résoud un système d’équations linéaires. On note : .. math:: \begin{array}{l} \mathbb{E} X = \frac{1}{n}\sum_{i=1}^n X_i \text{ et } \mathbb{E} Y = \frac{1}{n}\sum_{i=1}^n Y_i \\ \mathbb{E}(X^2) = \frac{1}{n}\sum_{i=1}^n X_i^2 \text{ et } \mathbb{E}(XY) = \frac{1}{n}\sum_{i=1}^n X_i Y_i \end{array} Finalement : .. math:: \begin{array}{l} a^* = \frac{ \mathbb{E}(XY) - \mathbb{E} X \mathbb{E} Y}{\mathbb{E}(X^2) - (\mathbb{E} X)^2} \text{ et } b^* = \mathbb{E} Y - a^* \mathbb{E} X \end{array} Q1 ~~ On génère un nuage de points avec le code suivant : .. code:: ipython3 import random def generate_xy(n=100, a=0.5, b=1): res = [] for i in range(0, n): x = random.uniform(0, 10) res.append((x, x*a + b + random.gauss(0,1))) return res generate_xy(10) .. parsed-literal:: [(2.996890181837922, 2.8750295096923186), (4.264526460045277, 2.324063943726332), (4.718648422500299, 3.052469543647318), (2.442689562115705, 3.861870829036456), (0.13558433730903707, 0.5754835901808546), (5.59230695209076, 1.6209924216651825), (7.610357428256408, 3.3202733390571373), (8.678403330137792, 4.96766236219644), (9.427259745518597, 6.385862058140737), (9.273956381823456, 4.938275166261537)] Q2 ~~ Ecrire une fonction qui calcule :math:`\mathbb{E} X, \mathbb{E} Y, \mathbb{E}(XY), \mathbb{E}(X^2)`. Plusieurs étudiants m’ont demandé ce qu’était ``obs``. C’est simplement le résultat de la fonction précédente. .. code:: ipython3 def calcule_exyxyx2(obs): sx = 0 sy = 0 sxy = 0 sx2 = 0 for x, y in obs: sx += x sy += y sxy += x * y sx2 += x * x n = len(obs) return sx/n, sy/n, sxy/n, sx2/n obs = generate_xy(10) calcule_exyxyx2(obs) .. parsed-literal:: (5.523441805914873, 3.850511796328412, 25.88928454527569, 38.98854258182378) Q3 ~~ Calculer les grandeurs :math:`a^*`, :math:`b^*`. A priori, on doit retrouver quelque chose d’assez proche des valeurs choisies pour la première question : :math:`a=0.5`, :math:`b=1`. .. code:: ipython3 def calcule_ab(obs): sx, sy, sxy, sx2 = calcule_exyxyx2(obs) a = (sxy - sx * sx) / (sx2 - sx**2) b = sy - a * sx return a, b calcule_ab(obs) .. parsed-literal:: (-0.5446995618974346, 6.859128128176218) Q4 ~~ Compléter le programme. .. code:: ipython3 import random def generate_caty(n=100, a=0.5, b=1, cats=["rouge", "vert", "bleu"]): res = [] for i in range(0, n): x = random.randint(0,2) cat = cats[x] res.append((cat, 10*x**2*a + b + random.gauss(0,1))) return res generate_caty(10) .. parsed-literal:: [('vert', 5.132157444058703), ('vert', 6.088324149707968), ('rouge', 0.16315983779393228), ('rouge', 0.9717657424738734), ('rouge', 2.843197432779423), ('rouge', 0.7204386278807904), ('bleu', 21.89226869979884), ('rouge', -0.16605748011543708), ('rouge', -0.02903894820027486), ('rouge', 0.5787816483863786)] Q5 ~~ On voudrait quand même faire une régression de la variable :math:`Y` sur la variable catégorielle. On construit une fonction qui assigne un numéro quelconque mais distinct à chaque catégorie. La fonction retourne un dictionnaire avec les catégories comme clé et les numéros comme valeurs. .. code:: ipython3 def numero_cat(obs): mapping = {} for color, y in obs: if color not in mapping: mapping[color] = len(mapping) return mapping obs = generate_caty(100) numero_cat(obs) .. parsed-literal:: {'bleu': 0, 'rouge': 2, 'vert': 1} Q6 ~~ On construit la matrice :math:`M_{ic}` tel que : :math:`M_{ic}` vaut 1 si :math:`c` est le numéro de la catégorie :math:`X_i`, 0 sinon. .. code:: ipython3 import numpy def construit_M(obs): mapping = numero_cat(obs) M = numpy.zeros((len(obs), 3)) for i, (color, y) in enumerate(obs): cat = mapping[color] M[i, cat] = 1.0 return M M = construit_M(obs) M[:5] .. parsed-literal:: array([[ 1., 0., 0.], [ 0., 1., 0.], [ 1., 0., 0.], [ 1., 0., 0.], [ 0., 0., 1.]]) Q7 ~~ Il est conseillé de convertir la matrice :math:`M` et les :math:`Y` au format *numpy*. On ajoute un vecteur constant à la matrice :math:`M`. La requête *numpy add column* sur un moteur de recherche vous directement à ce résultat : `How to add an extra column to an numpy array `__. .. code:: ipython3 def convert_numpy(obs): M = construit_M(obs) Mc = numpy.hstack([M, numpy.ones((M.shape[0], 1))]) Y = numpy.array([y for c, y in obs]) return M, Mc, Y.reshape((M.shape[0], 1)) M, Mc, Y = convert_numpy(obs) Mc[:5], Y[:5] .. parsed-literal:: (array([[ 1., 0., 0., 1.], [ 0., 1., 0., 1.], [ 1., 0., 0., 1.], [ 1., 0., 0., 1.], [ 0., 0., 1., 1.]]), array([[ 21.15485572], [ 6.37882494], [ 21.37124634], [ 21.77476221], [ 2.03305199]])) Q8 ~~ On résoud la régression multidimensionnelle en appliquant la formule :math:`C^* = (M'M)^{-1}M'Y`. La question 7 ne servait pas à grand chose excepté faire découvrir la fonction `hstack `__ car le rang de la matrice ``Mc`` est 3 < 4. .. code:: ipython3 alpha = numpy.linalg.inv(M.T @ M) @ M.T @ Y alpha .. parsed-literal:: array([[ 20.92499253], [ 6.14818418], [ 1.09988478]]) Q9 ~~ La régression détermine les coefficients :math:`\alpha` dans la régression :math:`Y_i = \alpha_{rouge} \mathbf{1\!\!1}_{X_i = rouge} + \alpha_{vert} \mathbf{1\!\!1}_{X_i = vert} + \alpha_{bleu} \mathbf{1\!\!1}_{X_i = bleu} + \epsilon_i`. Construire le vecteur :math:`\hat{Y_i} = \alpha_{rouge} \mathbf{1\!\!1}_{X_i = rouge} + \alpha_{vert} \mathbf{1\!\!1}_{X_i = vert} + \alpha_{bleu} \mathbf{1\!\!1}_{X_i = bleu}`. .. code:: ipython3 Yp = numpy.zeros((M.shape[0], 1)) for i in range(3): Yp[ M[:,i] == 1, 0] = alpha[i, 0] Yp[:5] .. parsed-literal:: array([[ 20.92499253], [ 6.14818418], [ 20.92499253], [ 20.92499253], [ 1.09988478]]) Q10 ~~~ Utiliser le résultat de la question 3 pour calculer les coefficients de la régression :math:`Y_i = a^* \hat{Y_i} + b^*`. .. code:: ipython3 obs = [(x, y) for x, y in zip(Yp, Y)] calcule_ab( obs ) .. parsed-literal:: (array([ 1.]), array([ 1.77635684e-15])) On aboutit au résultat :math:`Y = \hat{Y} + \epsilon`. On a associé une valeur à chaque catégorie de telle sorte que la régression de :math:`Y` sur cette valeur soit cette valeur. Autrement dit, c’est la meilleur approximation de :math:`Y` sur chaque catégorie. A quoi cela correspond-il ? C’est le second énoncé qui répond à cette question. Exercice 2 ---------- Q1 - Q2 - Q3 ~~~~~~~~~~~~ Ce sont les mêmes réponses. Q4 ~~ Un moyen très simple de simuler une `loi multinomiale `__ est de partir d’une loi uniforme et discrète à valeur dans entre 1 et 10. On tire un nombre, s’il est inférieur ou égal à 5, ce sera la catégorie 0, 1 si c’est inférieur à 8, 2 sinon. .. code:: ipython3 def generate_caty(n=100, a=0.5, b=1, cats=["rouge", "vert", "bleu"]): res = [] for i in range(0, n): # on veut 50% de rouge, 30% de vert, 20% de bleu x = random.randint(1, 10) if x <= 5: x = 0 elif x <= 8: x = 1 else : x = 2 cat = cats[x] res.append((cat, 10*x**2*a + b + random.gauss(0,1))) return res obs = generate_caty(10) obs .. parsed-literal:: [('vert', 6.0084452843428675), ('rouge', 2.155449750270483), ('rouge', 2.1132607428792447), ('vert', 6.897729973062269), ('rouge', 0.7637316114791164), ('vert', 5.566787193134299), ('vert', 5.848567708215508), ('bleu', 19.722503065860707), ('rouge', 0.8043492141543047), ('bleu', 21.675781652825997)] Q5 ~~ On voudrait quand même faire une régression de la variable :math:`Y` sur la variable catégorielle. On commence par les compter. Construire une fonction qui compte le nombre de fois qu’une catégorie est présente dans les données (un histogramme). .. code:: ipython3 def histogram_cat(obs): h = dict() for color, y in obs: h[color] = h.get(color, 0) + 1 return h histogram_cat(obs) .. parsed-literal:: {'bleu': 2, 'rouge': 4, 'vert': 4} Q6 ~~ Construire une fonction qui calcule la moyenne des :math:`Y_i` pour chaque catégorie : :math:`\mathbb{E}(Y | rouge)`, :math:`\mathbb{E}(Y | vert)`, :math:`\mathbb{E}(Y | bleu)`. La fonction retourne un dictionnaire ``{couleur:moyenne}``. .. code:: ipython3 def moyenne_cat(obs): h = dict() sy = dict() for color, y in obs: h[color] = h.get(color, 0) + 1 sy[color] = sy.get(color, 0) + y for k, v in h.items(): sy[k] /= v return sy moyenne_cat(obs) .. parsed-literal:: {'bleu': 20.69914235934335, 'rouge': 1.4591978296957873, 'vert': 6.080382539688736} L’énoncé induisait quelque peu en erreur car la fonction suggérée ne permet de calculer ces moyennes. Il suffit de changer. Q7 ~~ Construire le vecteur :math:`Z_i = \mathbb{E}(Y | rouge)\mathbf{1\!\!1}_{X_i = rouge} + \mathbb{E}(Y | vert) \mathbf{1\!\!1}_{X_i = vert} + \mathbb{E}(Y | bleu) \mathbf{1\!\!1}_{X_i = bleu}`. .. code:: ipython3 moys = moyenne_cat(obs) Z = [moys[c] for c, y in obs] Z[:5] .. parsed-literal:: [6.080382539688736, 1.4591978296957873, 1.4591978296957873, 6.080382539688736, 1.4591978296957873] Q8 ~~ Utiliser le résultat de la question 3 pour calculer les coefficients de la régression :math:`Y_i = a^* Z_i + b^*`. .. code:: ipython3 obs2 = [(z, y) for (c, y), z in zip(obs, Z)] calcule_ab( obs2 ) .. parsed-literal:: (1.0, 1.7763568394002505e-15) On aboutit au résultat :math:`Y = \hat{Y} + \epsilon`. On a associé une valeur à chaque catégorie de telle sorte que la régression de :math:`Y` sur cette valeur soit cette valeur. Q9 ~~ Calculer la matrice de variance / covariance pour les variables :math:`(Y_i)`, :math:`(Z_i)`, :math:`(Y_i - Z_i)`, :math:`\mathbf{1\!\!1}_{X_i = rouge}`, :math:`\mathbf{1\!\!1}_{X_i = vert}`, :math:`\mathbf{1\!\!1}_{X_i = bleu}`. .. code:: ipython3 bigM = numpy.empty((len(obs), 6)) bigM[:, 0] = [o[1] for o in obs] bigM[:, 1] = Z bigM[:, 2] = bigM[:, 0] - bigM[:, 1] bigM[:, 3] = [ 1 if o[0] == "rouge" else 0 for o in obs] bigM[:, 4] = [ 1 if o[0] == "vert" else 0 for o in obs] bigM[:, 5] = [ 1 if o[0] == "bleu" else 0 for o in obs] bigM[:5] .. parsed-literal:: array([[ 6.00844528, 6.08038254, -0.07193726, 0. , 1. , 0. ], [ 2.15544975, 1.45919783, 0.69625192, 1. , 0. , 0. ], [ 2.11326074, 1.45919783, 0.65406291, 1. , 0. , 0. ], [ 6.89772997, 6.08038254, 0.81734743, 0. , 1. , 0. ], [ 0.76373161, 1.45919783, -0.69546622, 1. , 0. , 0. ]]) On utilise la fonction `cov `__. .. code:: ipython3 c = numpy.cov(bigM.T) c .. parsed-literal:: array([[ 5.62221004e+01, 5.56972711e+01, 5.24829301e-01, -2.53176124e+00, -4.77901369e-01, 3.00966261e+00], [ 5.56972711e+01, 5.56972711e+01, -1.92890933e-16, -2.53176124e+00, -4.77901369e-01, 3.00966261e+00], [ 5.24829301e-01, -1.92890933e-16, 5.24829301e-01, -5.54535166e-17, 7.40725030e-17, -1.24510807e-17], [ -2.53176124e+00, -2.53176124e+00, -5.54535166e-17, 2.66666667e-01, -1.77777778e-01, -8.88888889e-02], [ -4.77901369e-01, -4.77901369e-01, 7.40725030e-17, -1.77777778e-01, 2.66666667e-01, -8.88888889e-02], [ 3.00966261e+00, 3.00966261e+00, -1.24510807e-17, -8.88888889e-02, -8.88888889e-02, 1.77777778e-01]]) On affiche un peu mieux les résultats : .. code:: ipython3 import pandas pandas.DataFrame(c).applymap(lambda x: '%1.3f' % x) .. raw:: html
0 1 2 3 4 5
0 56.222 55.697 0.525 -2.532 -0.478 3.010
1 55.697 55.697 -0.000 -2.532 -0.478 3.010
2 0.525 -0.000 0.525 -0.000 0.000 -0.000
3 -2.532 -2.532 -0.000 0.267 -0.178 -0.089
4 -0.478 -0.478 0.000 -0.178 0.267 -0.089
5 3.010 3.010 -0.000 -0.089 -0.089 0.178
Q10 ~~~ On permute rouge et vert. Construire le vecteur :math:`W_i = \mathbb{E}(Y | rouge)\mathbf{1\!\!1}_{X_i = vert} + \mathbb{E}(Y | vert)\mathbf{1\!\!1}_{X_i = rouge} + \mathbb{E}(Y | bleu)\mathbf{1\!\!1}_{X_i = bleu}`. Utiliser le résultat de la question 3 pour calculer les coefficients de la régression :math:`Y_i = a^* W_i + b^*`. Vérifiez que l’erreur est supérieure. .. code:: ipython3 moys = moyenne_cat(obs) moys["rouge"], moys["vert"] = moys.get("vert", 0), moys.get("rouge", 0) .. code:: ipython3 W = [moys[c] for c, y in obs] obs3 = [(w, y) for (c, y), w in zip(obs, W)] calcule_ab( obs3 ) .. parsed-literal:: (0.829591905722086, 1.2193824894893863) .. code:: ipython3 def calcule_erreur(obs): a, b = calcule_ab(obs) e = [(a*x + b - y)**2 for x, y in obs] return sum(e) / len(obs) calcule_erreur(obs2), calcule_erreur(obs3) .. parsed-literal:: (0.4723463712054069, 16.100975199731273) C’est carrément supérieur. Conclusion ---------- L’\ `analyse des correspondances multiples `__ est une façon d’étudier les modalités de variables catégorielles mais cela ne fait pas de la prédiction. Le modèle `logit - probit `__ prédit une variable binaire à partir de variables continue mais dans notre cas, c’est la variable à prédire qui est continue. Pour effectuer une prédiction, il convertit les catégories en variables numériques (voir `Categorical Variables `__). Le langage R est plus outillé pour cela : `Regression on categorical variables `__. Le module `categorical-encoding `__ est disponible en python. Cet examen décrit une méthode parmi d’autres pour transformer les catégories en variables continues.