Problème

Commencez par télécharger sur la page web du cours les jeux de données sur la page web du cours : http://www.di.ens.fr/appstat/TP/classification_data.zip

Les fichiers sont constitués de trois colonnes. Les deux premières donnent les coordonnées de points d’entrée $x_i ∈ ℝ$.

La troisième est la donnée de sortie correspondante $y_i ∈ \{0, 1\}$. A partir des données d’entraînement des fichiers .train, on veut construire des classifieurs qui seront efficaces sur les données des fichiers .test.

Génération des données

Le premier jeu de données est issu d’un modèle où les points de chaque classe sont générés suivant des loi Gaussiennes de variance commune entre classe mais de moyennes différentes, i.e, $P(X\mid Y = i)$ est une distribution Gaussienne de moyenne $µ_i$ et de matrice de covariance commune Σ.

Pour le deuxième jeu de données, le modèle de génération est le même mais cette fois les matrices de covariances aussi sont différentes. Concernant, le troisième jeu de données, une des classes est le mélange de deux Gaussiennes tandis que l’autre correspond à une seule Gaussienne.

Modèles de régression.

In [1]:
# /!\ Dans ce TP, j'utilise Python3 

import plotly
import numpy as np
import plotly.plotly as py
import plotly.graph_objs as go

import matplotlib.pyplot as plt
import scipy.io as sio
from scipy import stats
%matplotlib inline

plotly.offline.init_notebook_mode(connected=True)
from IPython.core.display import display, HTML
# The polling here is to ensure that plotly.js has already been loaded before
# setting display alignment in order to avoid a race condition.
display(HTML(
    '<script>'
        'var waitForPlotly = setInterval( function() {'
            'if( typeof(window.Plotly) !== "undefined" ){'
                'MathJax.Hub.Config({ SVG: { font: "STIX-Web" }, displayAlign: "center" });'
                'MathJax.Hub.Queue(["setRenderer", MathJax.Hub, "SVG"]);'
                'clearInterval(waitForPlotly);'
            '}}, 250 );'
    '</script>'
))

Commandes MathJax

$$ \newcommand{\e}{\mathop{\rm e}\nolimits} \newcommand{\Tr}{\mathop{\rm Tr}\nolimits} \newcommand{\GL}{\mathop{\rm GL}\nolimits} \newcommand{\Ker}{\mathop{\rm Ker}\nolimits} \newcommand{\I}{\mathop{\rm I}\nolimits} \newcommand{\Nor}{\mathop{\rm Nor}\nolimits} \newcommand{\Cen}{\mathop{\rm C}\nolimits} \newcommand{\Gr}{\mathop{\rm Gr}\nolimits} \newcommand{\ll}{\left[\!\!\left[} \newcommand{\rr}{\right]\!\!\right]} $$

1) Pour chacun des ensembles d’apprentissages A, B et C, représentez les données sous forme d’un nuage bicolore de points en 2D.

In [2]:
letters = ['A', 'B', 'C']
ranges = {}
data = {}
trace = {}
layout = {}

for letter in letters:
    
    data[letter] = {}
    trace[letter] = {}
    ranges[letter] = {}
    
    # Récupération des données d'entraînement et de test sous forme d'array
    for train_test in ['train', 'test']:
        
        with open('classification{}{}'.format(letter, train_test)) as f:
            lines = []
            
            for line in f:
                lines.append(list(map(float, line.split())))
            
            lines = np.array(lines).reshape([-1, 3])

            X, Y = lines[:, :2], lines[:, 2]
            
            data[letter][train_test] = {}
            data[letter][train_test]['X'] = X
            data[letter][train_test]['Y'] = Y
    
    X_1 = np.concatenate((data[letter]['train']['X'][:,0], data[letter]['test']['X'][:,0]))
    X_2 = np.concatenate((data[letter]['train']['X'][:,1], data[letter]['test']['X'][:,1]))
    
    ranges[letter]['x_1'] = [int(np.min(X_1)-1), int(np.max(X_1)+1)]
    ranges[letter]['x_2'] = [int(np.min(X_2)-2), int(np.max(X_2)+2)]
    
    trace[letter]['data'] = go.Scatter(
        x = X_1,
        y = X_2,
        mode = 'markers',
        marker = dict(
            size = 10,
            color = np.concatenate((
                2*data[letter]['train']['Y'], 
                2*data[letter]['test']['Y']+1
            )),
            colorbar = dict(
                title = 'y',
                titleside = 'top',
                tickmode = 'array',
                tickvals = [0, 1, 2, 3],
                ticktext = ['0 - entraînement','0 - test', '1 - entraînement', '1 - test'],
                ticks = 'outside'
            ),
            showscale = True,
            line = dict(
                width = 2,
                color = 'rgb(0, 0, 0)'
            ),
            colorscale = [
                [0, 'rgb(145,191,219)'], [0.25, 'rgb(145,191,219)'],
                [0.25, 'rgb(69,117,180)'], [0.5, 'rgb(69,117,180)'],
                [0.5, 'rgb(252,141,89)'], [0.75, 'rgb(252,141,89)'],
                [0.75, 'rgb(215,48,39)'], [1, 'rgb(215,48,39)']
            ]
        ),
        showlegend = False
    )

    layout[letter] = dict(
        title = 'Ensemble d\'apprentissage {}'.format(letter),
        xaxis = dict(
            title = '$x_1$',
            range = ranges[letter]['x_1'],
            ticklen = 5,
            zeroline = False,
            gridwidth = 2,
            ),
        yaxis = dict(
            title = '$x_2$',
            range = ranges[letter]['x_2'],
            ticklen = 5,
            gridwidth = 2,
            ),
        legend = dict(
            orientation = 'h',
            y = -0.2
        )
    )

    fig = dict(data=[trace[letter]['data']], layout=layout[letter])
    plotly.offline.iplot(fig)

2) Calculez le classifieur obtenu avec la régression linéaire.

In [3]:
classifiers = {}

classifiers['linear'] = {}

# Régression linéaire

def linear_regression(X, Y):
    # matrice de design 
    X = np.hstack((X, np.ones(len(X)).reshape([-1, 1])))

    # Estimateur des moindres carrés ordinaires (OLS): 
    # W = (X^T X)^{-1} X^T Y
    return np.linalg.solve(X.T.dot(X), X.T.dot(Y))

for letter in letters:

    X, Y = data[letter]['train']['X'], data[letter]['train']['Y']
    
    W = linear_regression(X, Y)
    
    # Stocker le classifier pour plus tard
    classifiers['linear'][letter] = W
    
    # Nombre d'erreurs sur l'ensemble de test
    
    a, b, c = W

    # Prédiction de la régression linéaire
    
    range_values = np.linspace(ranges[letter]['x_1'][0], ranges[letter]['x_1'][1], 20)

    trace[letter]['linear'] = go.Scatter(
        x = range_values,
        # a x_1 + b x_2 + c = 0.5 -> x_2 = (-a x_1 - c + 0.5)/b
        # ici: x = x_1, y = x_2 
        y = (-a*range_values - c + 0.5)/b,
        mode = 'lines',
        line = dict(
            color = ('rgb(7,142,145)'),
            width = 4,
            dash = 'dot'),
        name = 'Régression linéaire'
    )
    
    fig = dict(data=[trace[letter]['data'], trace[letter]['linear']], layout=layout[letter])
    plotly.offline.iplot(fig)

3) Implémentez le classifieur obtenu par régression logistique en utilisant soit une méthode de descente de gradient, soit l’algorithme des moindres carrés repondérés.

In [4]:
classifiers['logistic'] = {}

# cf. le cours http://www.di.ens.fr/appstat/notes/cours2-regression.pdf (pages 10-12)

def logistic_loss(y_hat, y):
    return y * np.log(1 + np.exp(-y_hat)) + (1 - y) * np.log(1 + np.exp(y_hat))

logistic_loss = np.vectorize(logistic_loss)

def sigmoid(x):
    return 1/(1 + np.exp(-x))

sigmoid = np.vectorize(sigmoid)

def risk(W, X, Y):
    return np.mean(logistic_loss(X.dot(W), Y))

def gradient_risk(W, X, Y, n):
    return (X.T.dot(sigmoid(X.dot(W)) - Y))/n

    
def logistic_regression(X, Y, display_risk=False):
    # matrice de design 
    X = np.hstack((X, np.ones(len(X)).reshape([-1, 1])))
    n = len(X)

    ################## Descente de gradient: début

    # Initialisation
    W = np.zeros(3)
    n_iter_max = 300
    eta = 10 # paramètre d'apprentissage
    eps = 10**(-4)  # critère d'arrêt de l'algorithme
    risks = [] # liste des risques calculés

    for i in range(n_iter_max):

        grad_W = gradient_risk(W, X, Y, n)

        # Mise à jour de W: W_{i+1} = W_i - eta * grad_W
        W -= eta * grad_W

        # Mise à jour du risque
        risks.append(risk(W, X, Y))

        if np.linalg.norm(grad_W) < eps:
            break

    if display_risk:
        # Affichage du risque lors de la descente de gradient 
        trace_risk = go.Scatter(
            x = list(range(1, len(risks)+1)),
            y = risks,
            mode = 'lines',
            name = 'Risque logistique'
        )

        layout_risk = go.Layout(
            title= 'Ensemble {} : Évolution du risque lors de la descente de gradient'.format(letter),
            hovermode= 'closest',
            xaxis= dict(
                title= 'Itérations',
                ticklen= 5,
                zeroline= False,
                gridwidth= 2,
            ),
            yaxis=dict(
                title= 'Risque',
                type='log',
                ticklen= 5,
                gridwidth= 2,
            ),
            showlegend= True
        )

        plotly.offline.iplot(go.Figure(data=[trace_risk], layout=layout_risk))

    ############# Descente de gradient: fin

    return W

for letter in letters:

    X, Y = data[letter]['train']['X'], data[letter]['train']['Y']

    W = logistic_regression(X, Y, display_risk=True)
    
    # Stocker le classifieur pour plus tard
    classifiers['logistic'][letter] = W

    a, b, c = W

    # Prédiction de la régression logistique
    
    range_values = np.linspace(ranges[letter]['x_1'][0], ranges[letter]['x_1'][1], 20)

    trace[letter]['logistic'] = go.Scatter(
        x = range_values,
        # sigmoid(a x_1 + b x_2 + c) = 0.5 -> x_2 = (-a x_1 - c)/b
        # ici: x = x_1, y = x_2 
        y = (-a*range_values - c)/b,
        mode = 'lines',
        line = dict(
            color = ('rgb(236,215,4)'),
            width = 4,
            dash = 'dash'),
        name = 'Régression logistique'
    )
    
    fig = dict(data=[trace[letter]['data'], trace[letter]['linear'], trace[letter]['logistic']],\
               layout=layout[letter])
    plotly.offline.iplot(fig)

Discriminant de Fisher

On propose un autre modèle de classification dit génératif. On modélise les données de chaque classe comme un nuage Gaussien de distribution $𝒩(µ_1, Σ)$ pour la classe $1$ (i.e., $Y = 1$) et $𝒩(µ_0, Σ)$ pour la classe $0$ (i.e., $Y = 0$), où la matrice de covariance $Σ$ est supposée la même pour les deux classes.

On pose aussi $π_1 = ℙ(Y = 1)$, la probabilité pour un exemple d’avoir le label $1$. Ce modèle s’appelle le modèle de l’analyse discriminante linéaire (LDA en anglais).

NB: Comme $π$ apparaît dans l'expression des Gaussiennes, j'ai changé l'énoncé et ai posé $$π_1 ≝ ℙ(Y = 1)$$

4)

a) Quel est le modèle conditionnel de $ℙ(Y |X = x)$ induit par ce modèle génératif ?

Pour tout $k∈ \{0, 1\}$, par la règle de Bayes :

$$ \begin{align} P(Y=k \mid X=x) &= \frac{P(Y=k) \,P(X=x \mid Y=k)}{P(X=x)} \\ &= \frac{P(Y=k) \, P(X=x \mid Y=k)}{\displaystyle\sum_{i=0}^1 P(Y=i)\, P(X=x \mid Y=i)} \\ &= \frac{1}{1 + \displaystyle\frac{P(Y = 1-k)\, P(X=x \mid Y = 1-k)}{P(Y = k)\, P(X=x \mid Y = k)}} \\ &= \frac{1}{1 + \exp\left(-\ln\left(\displaystyle\frac{P(Y = k)\, P(X=x \mid Y = k)}{P(Y = 1-k)\, P(X=x \mid Y = 1-k)}\right)\right)} \end{align} $$

$$ \begin{align} \ln\left(\displaystyle\frac{P(Y = k)\, P(X=x \mid Y = k)}{P(Y = 1-k)\, P(X=x \mid Y = 1-k)}\right) &= \ln\left(\displaystyle\frac{P(Y = k)}{P(Y = 1-k)}\right) + \ln\left(\frac{ \frac{1}{2π \sqrt{\det Σ}} \,\exp \left(- \frac 1 2 (x - μ_k)^T Σ^{-1} (x - μ_k) \right) }{ \frac{1}{2π \sqrt{\det Σ}} \,\exp \left(- \frac 1 2 (x - μ_{1-k})^T Σ^{-1} (x - μ_{1-k}) \right) }\right)\\ &= \ln\left(\displaystyle\frac{π_1^k (1-π_1)^{1-k}}{π_1^{1-k} (1-π_1)^k}\right) + \ln\left(\exp \left(-\frac 1 2 \Big((x - μ_k)^T Σ^{-1} (x - μ_k) - (x - μ_{1-k})^T Σ^{-1} (x - μ_{1-k})\Big) \right)\right) \\ &= \ln\left(\displaystyle\frac{π_1^k (1-π_1)^{1-k}}{π_1^{1-k} (1-π_1)^k}\right) -\frac 1 2 \Big( x^T Σ^{-1} x + μ_k^T Σ^{-1} μ_k \underbrace{- x^T Σ^{-1} μ_k - μ_k^T Σ^{-1} x}_{= \; - 2 x^T Σ^{-1} μ_k} - (x - μ_{1-k})^T Σ^{-1} (x - μ_{1-k})\Big) \\ &= \ln\left(\displaystyle\frac{π_1^k (1-π_1)^{1-k}}{π_1^{1-k} (1-π_1)^k}\right) -\frac 1 2 \Big(μ_k^T Σ^{-1} μ_k - μ_{1-k}^T Σ^{-1} μ_{1-k} - 2 x^T Σ^{-1} (μ_k - μ_{1-k}) \Big) \\ &= \underbrace{\ln\left(\displaystyle\frac{π_1^k (1-π_1)^{1-k}}{π_1^{1-k} (1-π_1)^k}\right) + \frac 1 2 \big(μ_{1-k}^T Σ^{-1} μ_{1-k} - μ_{k}^T Σ^{-1} μ_{k}\big)}_{≝ \; α_k} + x^T Σ^{-1} (μ_k - μ_{1-k}) \\ \end{align} $$

Donc

$$P(Y=k \mid X=x) = \frac{1}{1 + \e^{-\left(α_k + x^T Σ^{-1} (μ_k - μ_{1-k})\right)}}$$

où $$α_k ≝ \ln\left(\displaystyle\frac{π_1^k (1-π_1)^{1-k}}{π_1^{1-k} (1-π_1)^k}\right) + \frac 1 2 \big(μ_{1-k}^T Σ^{-1} μ_{1-k} - μ_{k}^T Σ^{-1} μ_{k}\big)$$

b) Montrez qu’il correspond à un modèle probabiliste de classification linéaire de la même forme que la régression logistique.

En posant $\tilde{x} = \begin{pmatrix}x \\ 1 \\ \end{pmatrix}$, il vient, avec $k=1$, que

$$P(Y=1 \mid X=x) = \frac{1}{1 + \exp\Big(-\big(\underbrace{α_1 + x^T Σ^{-1} (μ_1 - μ_0)}_{≝ \; f(\tilde{x})}\big)\Big)} = σ(f(\tilde{x}))$$

  • $$f(\tilde{x}) ≝ \tilde{x}^T \underbrace{\begin{pmatrix} Σ^{-1} (μ_1 - μ_0) \\ α_1 \\ \end{pmatrix}}_{≝ \; β}$$

  • $$α_1 ≝ \ln\left(\displaystyle\frac{π_1}{1-π_1}\right) + \frac 1 2 \big(μ_0^T Σ^{-1} μ_0 - μ_1^T Σ^{-1} μ_1\big)$$

$P(Y=1 \mid X=x)$ s'écrit donc sous la forme $σ(f(\tilde{x}))$, où $f(\tilde{x}) ≝ \tilde{x}^T β$ : l'analyse discriminante linéaire et la régression logistique utilisent le même modèle de classification linéaire sous-jacent.

c) Appliquez le modèle aux données en utilisant le principe du maximum de vraisemblance pour apprendre les paramètres du modèle génératif.

On procède de manière analogue à ce qu'on a vu dans le cours:

Principe du maximum de vraisemblance

  • $μ ≝ \begin{pmatrix} μ_0 \\ μ_1 \end{pmatrix} ∈ ℝ^2$
  • $Σ ∈ 𝒮_2(ℝ)$ est symmétrique définie positive

Dans la suite, on note

  • $x_i ≝ \begin{pmatrix} x_i^{(1)} \\ x_i^{(2)} \end{pmatrix}$ les points d'entrée
  • $\vec{x} ≝ (x_1, \ldots, x_n)$
  • $\vec{y} ≝ (y_1, \ldots, y_n)$

La log-vraissemblance s'écrit :

$$\begin{align*}

  • \ln L(\vec{x}, \vec{y}, μ, σ^2) &= -\sum\limits_{ i=1 }^n \ln \underbrace{P(X=x_i, Y=yi)}{=\; P(X=x_i \mid Y=y_i) P(Y=yi)}\ &= -\sum\limits{ i=1 }^n \ln \Big(\big(π_1 P(X=x_i \mid Y=1)\big)^{y_i}\big((1-π_1) P(X=x_i \mid Y=0)\big)^{1-yi}\Big)\ & = \underbrace{- \sum\limits{ i=1 }^n y_i \ln π_1 + (1-y_i) \ln(1-π_1) + y_i \ln\left(\frac{1}{2π \sqrt{\det Σ}} \exp \left(- \frac 1 2 (x_i - μ_1)^T Σ^{-1} (x_i - μ_1) \right)\right) + (1-y_i) \ln\left(\frac{1}{2π \sqrt{\det Σ}} \exp \left(- \frac 1 2 (x_i - μ_0)^T Σ^{-1} (x_i - μ0) \right)\right)}{≝ \; l(π_1, μ, Σ)} \ \end{align*}$$

Estimation de $π_1$

Vue comme une fonction de la variable $π_1$ seulement, $l$ est convexe et s'écrit sous la forme :

$$-\sum\limits_{ i=1 }^n y_i \ln π_1 + (1-y_i) \ln(1-π_1) + cste$$

L'annulation de la dérivée par rapport à $π_1$ donne :

$$ \begin{align*} & \sum\limits_{ i=1 }^n \frac{y_i}{π_1} - \frac{1-y_i}{1-π_1} = 0 \\ ⟺ \quad& π_1 = \frac{\sum\limits_{ i=1 }^n y_i}{\sum\limits_{ i=1 }^n y_i + (1-y_i)} \\ ⟺ \quad& π_1 = \frac 1 n \sum\limits_{ i=1 }^n y_i \end{align*} $$

$π_1$ est égal à la proportion de $x_i$ dont le label $y_i$ vaut $1$ (dans l'ensemble d'entraînement)

Estimation de $μ_1$

Vue comme une fonction de $μ_1$ seulement, $l$ est convexe et s'écrit sous la forme :

$$\sum\limits_{ i=1 }^n \frac{y_i}{2} (x_i - μ_1)^T Σ^{-1} (x_i - μ_1) + cste$$

L'annulation de la dérivée par rapport à $μ_1$ donne :

$$ \begin{align*} & \sum\limits_{ i=1 }^n y_i Σ^{-1} (x_i - μ_1) = 0 \\ ⟺ \quad& μ_1 = \frac{\sum\limits_{ i=1 }^n y_i x_i}{\sum\limits_{ i=1 }^n y_i} \end{align*} $$

De même :

$$μ_0 = \frac{\sum\limits_{ i=1 }^n (1-y_i) x_i}{\sum\limits_{ i=1 }^n (1-y_i)}$$

$μ_1$ (resp. $μ_0$) est égal à la moyenne des $x_i$ dont le label $y_i$ vaut $1$ (resp. $0$) (dans l'ensemble d'entraînement).

Estimation de $Σ$

Vue comme une fonction de $Λ ≝ Σ^{-1}$ seulement, $l$ est convexe et s'écrit sous la forme :

$$\sum\limits_{ i=1 }^n -\frac{1}{2} \ln(\vert Λ \vert) + \frac{y_i}{2} \underbrace{\underbrace{(x_i - μ_1)^T}_{1 × 2} \underbrace{Λ}_{2 × 2} \underbrace{(x_i - μ_1)}_{2 × 1}}_{\Tr((x_i - μ_1)^T Λ (x_i - μ_1))} + \frac{1-y_i}{2} \underbrace{\underbrace{(x_i - μ_0)^T}_{1 × 2} \underbrace{Λ}_{2 × 2} \underbrace{(x_i - μ_0)}_{2 × 1}}_{=\, \Tr((x_i - μ_0)^T Λ (x_i - μ_0))} + cste$$

L'annulation de la dérivée par rapport à $Λ$ donne :

$$ \begin{align*} & \sum\limits_{ i=1 }^n - Λ^{-1} + y_i(x_i - μ_1) (x_i - μ_1)^T + (1-y_i)(x_i - μ_0) (x_i - μ_0)^T = 0 \\ ⟺ \quad& Σ = Λ^{-1} = \frac 1 n \sum\limits_{ i=1 }^n y_i(x_i - μ_1) (x_i - μ_1)^T + (1-y_i)(x_i - μ_0) (x_i - μ_0)^T \end{align*} $$


Frontière de classification

On prédit $y=1$ si, et seulement si :

$$ \begin{align*} & P(X=x, Y=1) ≥ P(X=x, Y=0) \\ ⟺ \quad& \frac{P(X=x, Y=1)}{P(X=x, Y=0)} ≥ 1 \\ ⟺ \quad& \ln \left(\frac{P(X=x, Y=1)}{P(X=x, Y=0)}\right) ≥ 0 \\ ⟺ \quad& \ln \left(\frac{P(X=x \mid Y=1) P(Y=1)}{P(X=x \mid Y=0) P(Y=0)}\right) ≥ 0 \\ ⟺ \quad& \underbrace{\ln\left(\displaystyle\frac{π_1}{1-π_1}\right) + \frac 1 2 \big(μ_0^T Σ^{-1} μ_0 - μ_1^T Σ^{-1} μ_1\big)}_{≝ \; α_1} + x^T Σ^{-1} (μ_1 - μ_0) ≥ 0 &&\text{(cf. question précédente)}\\ ⟺ \quad& \tilde{x}^T \underbrace{\begin{pmatrix} Σ^{-1} (μ_1 - μ_0) \\ α_1 \\ \end{pmatrix}}_{≝ \; β} ≥ 0 &&\text{(avec les notations de la question précédente)} \end{align*} $$

In [5]:
classifiers['LDA'] = {}

# LDA : 
# classification binaire : les données sont catégorisées en 2 classes :
# 1. les x_i dont le label y_i vaut c
# 2. les autres x_i
# Ici : c=1 (l'autre classe étant 0)
# NB: fait de procéder ainsi facilitera les choses plus tard, pour MNIST

def LDA(X, Y, c=1):
    n = len(X)
    Y = Y.flatten()
    
    # pi_1 : proportion de x_i dont le label y_i vaut c
    pi_1 = np.mean(Y == c)
    
    # mu_1 (resp. mu_0) : vecteur moyen des x_i dont le label y_i vaut c (resp. non c) 
    mu_1 = np.mean(X[(Y == c), :], axis=0)
    mu_0 = np.mean(X[~(Y == c), :], axis=0)
    
    # Sigma : 
    Sigma = np.sum((x-mu_1).reshape([-1, 1]).dot((x-mu_1).reshape([1, -1])) 
                        if y == c else
                   (x-mu_0).reshape([-1, 1]).dot((x-mu_0).reshape([1, -1])) 
                        for x, y in zip(X, Y))/n
    
    # alpha_1, beta
    SigmaInv = np.linalg.pinv(Sigma) # Sigma^{-1}
    
    alpha_1 = np.log(pi_1/(1-pi_1)) \
            + mu_0.dot(SigmaInv.dot(mu_0))/2 - mu_1.dot(SigmaInv.dot(mu_1))/2
    
    beta = np.concatenate((SigmaInv.dot(mu_1 - mu_0), np.array([alpha_1])))
    
    return beta


for letter in letters:

    X, Y = data[letter]['train']['X'], data[letter]['train']['Y']
    
    beta = LDA(X, Y)
    
    # matrice de design 
    X = np.hstack((X, np.ones(len(X)).reshape([-1, 1])))
    
    # Stocker le classifieur pour plus tard
    classifiers['LDA'][letter] = beta

    a, b, c = beta

    # Prédiction de l'analyse discriminante linéaire
    
    range_values = np.linspace(ranges[letter]['x_1'][0], ranges[letter]['x_1'][1], 20)

    trace[letter]['LDA'] = go.Scatter(
        x = range_values,
        # a x_1 + b x_2 + c = 0 -> x_2 = (-a x_1 - c)/b
        # ici: x = x_1, y = x_2 
        y = (-a*range_values - c)/b,
        mode = 'lines',
        line = dict(
            color = ('rgb(221,52,151)'),
            width = 1),
        name = 'LDA'
    )
    
    fig = dict(data=[trace[letter]['data'], trace[letter]['linear'], trace[letter]['logistic'],\
                     trace[letter]['LDA']], layout=layout[letter])
    plotly.offline.iplot(fig)