1A.e - TD noté, 16 décembre 2016#

Régression linéaire avec des variables catégorielles.

Exercice 1#

On suppose qu’on dispose d’un ensemble d’observations $(X_i, Y_i)$ avec $X_i, Y_i \in \mathbb{R}$ . La régression linéaire consiste une relation linéaire $Y_i = a X_i + b + \epsilon_i$ qui minimise la variance du bruit. On pose :

$E(a, b) = \sum_i (Y_i - (a X_i + b))^2$

On cherche $a, b$ tels que :

$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 :

$\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 :

$\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 :

$\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 :

[2]:

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)

[2]:

[(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 $\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.

[3]:

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)

[3]:

(5.523441805914873, 3.850511796328412, 25.88928454527569, 38.98854258182378)

Q3#

Calculer les grandeurs $a^*$ , $b^*$ . A priori, on doit retrouver quelque chose d’assez proche des valeurs choisies pour la première question : $a=0.5$ , $b=1$ .

[4]:

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)

[4]:

(-0.5446995618974346, 6.859128128176218)

Q4#

Compléter le programme.

[5]:

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)

[5]:

[('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 $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.

[6]:

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)

[6]:

{'bleu': 0, 'rouge': 2, 'vert': 1}

Q6#

On construit la matrice $M_{ic}$ tel que : $M_{ic}$ vaut 1 si $c$ est le numéro de la catégorie $X_i$ , 0 sinon.

[7]:

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]

[7]:

array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 1.,  0.,  0.],
       [ 1.,  0.,  0.],
       [ 0.,  0.,  1.]])

Q7#

Il est conseillé de convertir la matrice $M$ et les $Y$ au format numpy. On ajoute un vecteur constant à la matrice $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.

[8]:

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]

[8]:

(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 $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.

[9]:

alpha = numpy.linalg.inv(M.T @ M) @ M.T @ Y
alpha

[9]:

array([[ 20.92499253],
       [  6.14818418],
       [  1.09988478]])

Q9#

La régression détermine les coefficients $\alpha$ dans la régression $Y_i = \alpha_{rouge} \mathbb{1}_{X_i = rouge} + \alpha_{vert} \mathbb{1}_{X_i = vert} + \alpha_{bleu} \mathbb{1}_{X_i = bleu} + \epsilon_i$ .

Construire le vecteur $\hat{Y_i} = \alpha_{rouge} \mathbb{1}_{X_i = rouge} + \alpha_{vert} \mathbb{1}_{X_i = vert} + \alpha_{bleu} \mathbb{1}_{X_i = bleu}$ .

[10]:

Yp = numpy.zeros((M.shape[0], 1))
for i in range(3):
    Yp[M[:, i] == 1, 0] = alpha[i, 0]
Yp[:5]

[10]:

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 $Y_i = a^* \hat{Y_i} + b^*$ .

[11]:

obs = [(x, y) for x, y in zip(Yp, Y)]
calcule_ab(obs)

[11]:

(array([ 1.]), array([  1.77635684e-15]))

On aboutit au résultat $Y = \hat{Y} + \epsilon$ . On a associé une valeur à chaque catégorie de telle sorte que la régression de $Y$ sur cette valeur soit cette valeur. Autrement dit, c’est la meilleur approximation de $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.

[12]:

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

[12]:

[('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 $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).

[13]:

def histogram_cat(obs):
    h = dict()
    for color, y in obs:
        h[color] = h.get(color, 0) + 1
    return h


histogram_cat(obs)

[13]:

{'bleu': 2, 'rouge': 4, 'vert': 4}

Q6#

Construire une fonction qui calcule la moyenne des $Y_i$ pour chaque catégorie : $\mathbb{E}(Y | rouge)$ , $\mathbb{E}(Y | vert)$ , $\mathbb{E}(Y | bleu)$ . La fonction retourne un dictionnaire {couleur:moyenne}.

[14]:

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)

[14]:

{'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 $Z_i = \mathbb{E}(Y | rouge)\mathbb{1}_{X_i = rouge} + \mathbb{E}(Y | vert) \mathbb{1}_{X_i = vert} + \mathbb{E}(Y | bleu) \mathbb{1}_{X_i = bleu}$ .

[15]:

moys = moyenne_cat(obs)
Z = [moys[c] for c, y in obs]
Z[:5]

[15]:

[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 $Y_i = a^* Z_i + b^*$ .

[16]:

obs2 = [(z, y) for (c, y), z in zip(obs, Z)]
calcule_ab(obs2)

[16]:

(1.0, 1.7763568394002505e-15)

On aboutit au résultat $Y = \hat{Y} + \epsilon$ . On a associé une valeur à chaque catégorie de telle sorte que la régression de $Y$ sur cette valeur soit cette valeur.

Q9#

Calculer la matrice de variance / covariance pour les variables $(Y_i)$ , $(Z_i)$ , $(Y_i - Z_i)$ , $\mathbb{1}_{X_i = rouge}$ , $\mathbb{1}_{X_i = vert}$ , $\mathbb{1}_{X_i = bleu}$ .

[17]:

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]

[17]:

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.

[18]:

c = numpy.cov(bigM.T)
c

[18]:

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 :

[19]:

import pandas

pandas.DataFrame(c).applymap(lambda x: "%1.3f" % x)

[19]:

	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 $W_i = \mathbb{E}(Y | rouge)\mathbb{1}_{X_i = vert} + \mathbb{E}(Y | vert)\mathbb{1}_{X_i = rouge} + \mathbb{E}(Y | bleu)\mathbb{1}_{X_i = bleu}$ . Utiliser le résultat de la question 3 pour calculer les coefficients de la régression $Y_i = a^* W_i + b^*$ . Vérifiez que l’erreur est supérieure.

[20]:

moys = moyenne_cat(obs)
moys["rouge"], moys["vert"] = moys.get("vert", 0), moys.get("rouge", 0)

[21]:

W = [moys[c] for c, y in obs]
obs3 = [(w, y) for (c, y), w in zip(obs, W)]
calcule_ab(obs3)

[21]:

(0.829591905722086, 1.2193824894893863)

[22]:

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)

[22]:

(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.

[23]:

Notebook on github