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 :

[1]:

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)

[1]:

[(6.049169429448925, 4.718935755862216),
 (4.627795365365467, 3.8048865083718906),
 (6.474192602737176, 5.038157546858375),
 (5.257002270747765, 3.824981234020261),
 (3.6776757481314792, 2.3846063467327263),
 (4.806391225454085, 3.6643444372898646),
 (8.524342428472508, 3.7743268407183237),
 (8.293296420000031, 6.157037850737982),
 (1.1448844698514993, 1.8337680880279141),
 (0.47468792840570706, 1.3511445685085155)]

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.

[2]:

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)

[2]:

(5.486540369984782, 3.9453189049350192, 25.815898347439713, 39.24240425818958)

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

[3]:

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)

[3]:

(-0.46893829734349396, 6.518167804342026)

Q4¶

Compléter le programme.

[4]:

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)

[4]:

[('bleu', 19.91792858279941),
 ('vert', 6.552945331178698),
 ('bleu', 21.607291021269408),
 ('vert', 6.299785531549073),
 ('vert', 7.258812927878141),
 ('rouge', -0.6578105813916351),
 ('rouge', 0.522136653538479),
 ('bleu', 21.469646109623795),
 ('rouge', 1.3477906240590773),
 ('vert', 7.922241071711617)]

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.

[5]:

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)

[5]:

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

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.

[6]:

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]

[6]:

array([[1., 0., 0.],
       [0., 1., 0.],
       [1., 0., 0.],
       [0., 0., 1.],
       [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.

[7]:

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]

[7]:

(array([[1., 0., 0., 1.],
        [0., 1., 0., 1.],
        [1., 0., 0., 1.],
        [0., 0., 1., 1.],
        [0., 0., 1., 1.]]),
 array([[ 4.7323734 ],
        [21.71040171],
        [ 4.84924589],
        [ 1.78526719],
        [ 0.34311133]]))

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.

[8]:

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

[8]:

array([[ 5.94180865],
       [20.97848259],
       [ 1.05427197]])

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

[9]:

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

[9]:

array([[ 5.94180865],
       [20.97848259],
       [ 5.94180865],
       [ 1.05427197],
       [ 1.05427197]])

Q10¶

Utiliser le résultat de la question 3 pour calculer les coefficients de la régression $Y_i = a^* \hat{Y_i} + b^*$ .

[10]:

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

[10]:

(array([1.]), array([-5.32907052e-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.

[11]:

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

[11]:

[('rouge', 0.09329495469271787),
 ('rouge', 0.789609852660202),
 ('vert', 4.714873982539039),
 ('rouge', 2.0829663032267485),
 ('vert', 6.745691356870579),
 ('rouge', 0.05932073340343336),
 ('vert', 6.772334898512818),
 ('vert', 5.8598815352891425),
 ('rouge', 2.30150566164451),
 ('vert', 6.549148011451264)]

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

[12]:

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


histogram_cat(obs)

[12]:

{'rouge': 5, 'vert': 5}

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

[13]:

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)

[13]:

{'rouge': 1.0653395011255224, 'vert': 6.1283859569325685}

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

[14]:

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

[14]:

[1.0653395011255224,
 1.0653395011255224,
 6.1283859569325685,
 1.0653395011255224,
 6.1283859569325685]

Q8¶

Utiliser le résultat de la question 3 pour calculer les coefficients de la régression $Y_i = a^* Z_i + b^*$ .

[15]:

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

[15]:

(1.0, -8.881784197001252e-16)

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

[16]:

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]

[16]:

array([[ 0.09329495,  1.0653395 , -0.97204455,  1.        ,  0.        ,
         0.        ],
       [ 0.78960985,  1.0653395 , -0.27572965,  1.        ,  0.        ,
         0.        ],
       [ 4.71487398,  6.12838596, -1.41351197,  0.        ,  1.        ,
         0.        ],
       [ 2.0829663 ,  1.0653395 ,  1.0176268 ,  1.        ,  0.        ,
         0.        ],
       [ 6.74569136,  6.12838596,  0.6173054 ,  0.        ,  1.        ,
         0.        ]])

On utilise la fonction cov.

[17]:

c = numpy.cov(bigM.T)
c

[17]:

array([[ 7.96951427e+00,  7.12067761e+00,  8.48836657e-01,
        -1.40640179e+00,  1.40640179e+00,  0.00000000e+00],
       [ 7.12067761e+00,  7.12067761e+00,  2.25880946e-16,
        -1.40640179e+00,  1.40640179e+00,  0.00000000e+00],
       [ 8.48836657e-01,  2.25880946e-16,  8.48836657e-01,
        -4.00913870e-17,  4.00913870e-17,  0.00000000e+00],
       [-1.40640179e+00, -1.40640179e+00, -4.00913870e-17,
         2.77777778e-01, -2.77777778e-01,  0.00000000e+00],
       [ 1.40640179e+00,  1.40640179e+00,  4.00913870e-17,
        -2.77777778e-01,  2.77777778e-01,  0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
         0.00000000e+00,  0.00000000e+00,  0.00000000e+00]])

On affiche un peu mieux les résultats :

[18]:

import pandas

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

[18]:

	0	1	2	3	4
0	7.970	7.121	0.849	-1.406	1.406
1	7.121	7.121	0.000	-1.406	1.406
2	0.849	0.000	0.849	-0.000	0.000
3	-1.406	-1.406	-0.000	0.278	-0.278
4	1.406	1.406	0.000	-0.278	0.278
5	0.000	0.000	0.000	0.000	0.000

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.

[19]:

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

[20]:

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

[20]:

(-0.999999999999999, 7.193725458058086)

[21]:

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)

[21]:

(0.7639529913942368, 0.7639529913942368)

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.

[ ]:

Notebook on github