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{split}\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}\end{split}\]

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