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.
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.
# /!\ 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>'
))
$$ \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]} $$
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)
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)
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)
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)$$
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} $$
où
$$ \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)$$
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}))$$
où
$$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.
On procède de manière analogue à ce qu'on a vu dans le cours:
Dans la suite, on note
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=y_i)}_{=\; P(X=x_i \mid Y=y_i) P(Y=y_i)}\\ &= -\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-y_i}\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*}$$
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)
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).
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*} $$
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*} $$
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)
Si on autorise les matrices de covariances à être différentes pour les deux classes (on les note $Σ_0$ et $Σ_1$), il vient, en reprenant les calculs, que :
On prédit $y=1$ si, et seulement si :
$$ \begin{align*} & \underbrace{\ln\left(\displaystyle\frac{π_1}{1-π_1}\right) + \frac 1 2 \big(μ_0^T Σ_0^{-1} μ_0 - μ_1^T Σ_1^{-1} μ_1\big)}_{≝ \; α'_1} + x^T (Σ_1^{-1} μ_1 - Σ_0^{-1} μ_0) + x^T (Σ_0^{-1} - Σ_1^{-1}) x ≥ 0 \\ ⟺ \quad& \tilde{x}^T \underbrace{\begin{pmatrix} Σ_1^{-1} μ_1 - Σ_0^{-1} μ_0 \\ α'_1 \\ \end{pmatrix}}_{≝ \; β'} + x^T (Σ_0^{-1} - Σ_1^{-1}) x ≥ 0 &&\text{(avec les notations des questions précédentes)} \end{align*} $$
D'où l'adjectif "quadratique" : la frontière de classification a un terme d'ordre $2$.
classifiers['QDA'] = {}
def QDA(X, Y, c=1):
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_k :
Sigma_1 = np.sum((x-mu_1).reshape([-1, 1]).dot((x-mu_1).reshape([1, -1]))\
for x, y in zip(X, Y) if y == c)/np.sum(Y == c)
Sigma_0 = np.sum((x-mu_0).reshape([-1, 1]).dot((x-mu_0).reshape([1, -1]))\
for x, y in zip(X, Y) if y != c)/np.sum(~(Y == c))
# alpha_prime_1, beta_prime
Sigma1Inv, Sigma0Inv = np.linalg.pinv(Sigma_1), np.linalg.pinv(Sigma_0) # Sigma_k^{-1}
alpha_prime_1 = np.log(pi_1/(1-pi_1)) \
+ mu_0.dot(Sigma0Inv.dot(mu_0))/2 - mu_1.dot(Sigma1Inv.dot(mu_1))/2
beta_prime = np.concatenate((Sigma1Inv.dot(mu_1) - Sigma0Inv.dot(mu_0), np.array([alpha_prime_1])))
# pour le terme d'ordre 2 : Sigma_0^{-1} - Sigma_1^{-1}
SigmaInvMinus = Sigma0Inv - Sigma1Inv
return beta_prime, SigmaInvMinus
for letter in letters:
X, Y = data[letter]['train']['X'], data[letter]['train']['Y']
beta_prime, SigmaInvMinus = QDA(X, Y)
# matrice de design
X = np.hstack((X, np.ones(len(X)).reshape([-1, 1])))
# Stocker le classifieur pour plus tard
classifiers['QDA'][letter] = beta_prime, SigmaInvMinus
a, b, c = beta_prime
d, e, f, g = SigmaInvMinus.flatten()
# Prédiction de l'analyse discriminante quadratique
x1_range_values = np.linspace(ranges[letter]['x_1'][0], ranges[letter]['x_1'][1], 100)
x2_range_values = np.linspace(ranges[letter]['x_2'][0], ranges[letter]['x_2'][1], 100)
x_1, x_2 = np.meshgrid(x1_range_values, x2_range_values)
# a x_1 + b x_2 + c + x^T (Sigma_0^{-1} - Sigma_1^{-1}) x = 0
contour_lines = a*x_1 + b*x_2 + c + x_1*(d*x_1+e*x_2) + x_2*(f*x_1+g*x_2)
trace[letter]['QDA'] = go.Contour(
x = x1_range_values,
y = x2_range_values,
z = contour_lines,
contours=dict(
start = 0,
end = 0,
),
colorbar = dict(
title = 'QDA',
x = -0.17,
tickmode = 'array',
tickvals = [0],
ticktext = ['y prédit'],
tickangle = -90,
ticks = 'inside',
),
colorscale = [[0, 'rgba(69,117,180,.5)'], [1, 'rgba(215,48,39,.5)']]
)
fig = dict(data=[trace[letter]['data'], trace[letter]['QDA']], layout=layout[letter])
plotly.offline.iplot(fig)
for letter in letters:
fig = dict(data=[trace[letter]['data'], trace[letter]['linear'], trace[letter]['logistic'],\
trace[letter]['LDA'], trace[letter]['QDA']], layout=layout[letter])
plotly.offline.iplot(fig)
from IPython.display import HTML
# Taux d'erreur
error_rate = {}
for letter in letters:
error_rate[letter] = {}
def compute_error_rates(letters, list_classifiers, error_rate=error_rate):
for letter in letters:
#### Pour chaque classifieur
for classifier in list_classifiers:
error_rate[letter][classifier] = {}
# Taux d'erreur sur l'ensemble d'entraînement et de test
for train_test in ['train', 'test']:
X, Y = data[letter][train_test]['X'], data[letter][train_test]['Y']
Y = Y.flatten()
# matrice de design
X = np.hstack((X, np.ones(len(X)).reshape([-1, 1])))
W = classifiers[classifier][letter]
# taux d'erreur
if classifier == 'linear':
error_rate[letter][classifier][train_test] = np.mean(((X.dot(W) >= 0.5) != Y))
elif classifier == 'logistic':
error_rate[letter][classifier][train_test] = np.mean(((sigmoid(X.dot(W)) >= 0.5) != Y))
elif classifier == 'LDA':
if letter != 'mnist':
error_rate[letter][classifier][train_test] = np.mean((X.dot(W) >= 0.) != Y)
else:
# Pour plus tard, avec MNIST :
# pour tout x_j :
# le W[i] qui maximise x_j W[i] détermine la classe (le chiffre) prédite de x_j
error_rate[letter][classifier][train_test] = 0
for ind, x in enumerate(X):
maximized_class = np.argmax(np.array([x.dot(w) for w in W]))
if maximized_class != Y[ind]:
error_rate[letter][classifier][train_test] += 1
error_rate[letter][classifier][train_test] /= len(X)
elif classifier == 'QDA':
if letter != 'mnist':
frontier = np.array([x.dot(W[0]) + x[:-1].dot(W[1].dot(x[:-1])) >= 0 for x in X])
error_rate[letter][classifier][train_test] = np.mean(frontier != Y)
else:
# Pour plus tard, avec MNIST :
error_rate[letter][classifier][train_test] = 0
for ind, x in enumerate(X):
frontier = np.array([x.dot(w[0]) + x[:-1].dot(w[1].dot(x[:-1])) for w in W])
maximized_class = np.argmax(frontier)
if maximized_class != Y[ind]:
error_rate[letter][classifier][train_test] += 1
error_rate[letter][classifier][train_test] /= len(X)
return error_rate
error_rate = compute_error_rates(letters, ['linear', 'logistic', 'LDA', 'QDA'])
# Affichage du tableau des taux d'erreurs
s = """<table>
<tr>
<th style='text-align:center'>A</th>
<th style='text-align:center'>B</th>
<th style='text-align:center'>C</th>
</tr>
<tr>
"""
for letter in letters:
s += """<td>
<table>
<tr>
<th></th>
<th style='text-align:center'>train</th>
<th style='text-align:center'>test</th>
</tr>
"""
for classifier in classifiers:
s += """<tr>
<th>{}</th>
<td>{}</td>
<td>{}</td>
""".format(classifier, error_rate[letter][classifier]['train'],\
error_rate[letter][classifier]['test'])
s += """</tr>
</table>
</td>
"""
s += """</tr>
</table>
"""
HTML(s)
# Diagramme en bâtons
trace_bar = {}
def bar_chart_error(letters, list_classifiers, list_colors, trace_bar = trace_bar):
learning_sets = ['{} {}'.format(letter, train_test)
for letter in letters for train_test in ['train', 'test']]
for classifier, color in zip(list_classifiers, list_colors):
trace_bar[classifier] = go.Bar(
x = learning_sets,
y = [error_rate[letter][classifier][train_test]
for letter in letters for train_test in ['train', 'test']],
name = classifier,
marker=dict(
color=color,
),
opacity=0.75
)
return trace_bar
list_classifiers = ['linear', 'logistic', 'LDA', 'QDA']
list_colors = ['rgb(7,142,145)', 'rgb(236,215,4)', 'rgb(221,52,151)', 'rgb(215,48,39)']
trace_bar = bar_chart_error(letters, list_classifiers, list_colors)
layout_bar = go.Layout(
barmode='group',
title = 'Taux d\'erreur : comparaison',
xaxis = dict(
title='Ensembles d\'apprentissage'
),
yaxis = dict(
title = 'Taux d\'erreur'
),
)
fig = go.Figure(data=[trace_bar[key] for key in trace_bar], layout=layout_bar)
plotly.offline.iplot(fig)
Dans les différents ensembles d'apprentissage, je parlerai de "points bleus" (resp. "points rouges") pour désigner les points dont le label est 0 (resp. 1).
Dans l'ensemble A : les points bleus et rouges semblent générés par deux Gaussiennes (une pour chaque couleur) dont la matrice de covariance a un valeur propre bien plus grande que l'autre, donnant une forme d'"ellispe allongée" aux amas bleus et rouge. De plus, ils se chevauchent un peu, mais semblent néanmoins séparables par une droite, d'où :
AOn notera tout de même que QDA est meilleur que les autres sur l'ensemble de test : la conique de séparation épouse un peu la frontière de séparation des amas bleus et rouges.
Dans l'ensemble B : les amas bleus et rouges semblent encore générés par deux Gaussiennes :
AIl est bien plus dur de séparer linéairement les deux amas, puisque l'amas rouge "plonge" (s'enfonce) dans l'amas bleu, d'où un taux d'erreur des classifieurs plus élevé que précédemment.
la régression logistique et QDA se démarque en étant meilleurs sur l'ensemble d'entraînement (la régression logisitique étant aussi meilleure que les autres sur l'ensemble test):
L'ensemble C est très similaire à B, à l'exception d'un deuxième amas circulaire de points rouges isolés de l'amas principal. Cet amas d'"outliers" fait baisser la performance de
Par contre, la frontière de la régression logistique passe en-dessous de l'amas : elle ne change donc pas, et reste plus performante que les autres.
from scipy.io import loadmat
mnist = loadmat('mnist_digits.mat') # données sous forme de dictionnaire
data['mnist'] = {}
######### Génération des données : début
nbr_img = 6000 # nombre d'images à sélectionner au hasard
ind_train = range(0, int(1/3 * nbr_img))
ind_test = range(int(1/3 * nbr_img), nbr_img)
ind_rand = np.random.randint(len(mnist['x']), size=nbr_img) # indices des images sélectionnées
data['mnist']['train'], data['mnist']['test'] = {}, {}
X_mnist, Y_mnist = mnist['x'][ind_rand, :], mnist['y'][ind_rand] # images, labels
data['mnist']['train']['X'], data['mnist']['train']['Y'] = X_mnist[ind_train, :], Y_mnist[ind_train]
data['mnist']['test']['X'], data['mnist']['test']['Y'] = X_mnist[ind_test, :], Y_mnist[ind_test]
# En considérant l'expression de la/les matrice(s) Sigma(_k) de la LDA (QDA),
# on se rend compte que : s'il existe un chiffre k tel qu'il existe une paire de
# pixels qui sont d'une couleur c fixée dans toutes les images représentant k,
# alors Sigma n'est pas inversible, puisqu'elle a deux colonnes liées
# Pour éviter ce cas de figure, on calcule le pseudo-inverse de
# la matrice de covariance, dans LDA/QDA
########### Génération des données : fin
#### Généralisation de LDA, QDA à une classification à nb_classes classes : 0, ..., nb_classes-1
## : on fait nb_classes classifications binaires : entre la classe 0 <= c < nb_classes et la classe non c
def LDA_generalized(X, Y, nb_classes = 10):
return [LDA(X, Y, c=c) for c in range(nb_classes)]
def QDA_generalized(X, Y, nb_classes = 10):
return [QDA(X, Y, c=c) for c in range(nb_classes)]
list_classifiers = ['LDA', 'QDA']
classifier_functions = {'LDA': LDA_generalized, 'QDA': QDA_generalized}
list_colors = ['rgb(221,52,151)', 'rgb(215,48,39)']
#### Calcul des classifieurs LDA et QDA pour MNIST
for classifier in list_classifiers:
classifiers[classifier]['mnist'] = classifier_functions[classifier](data['mnist']['train']['X'],\
data['mnist']['train']['Y'])
error_rate['mnist'] = {}
error_rate = compute_error_rates(letters+['mnist'], list_classifiers)
trace_bar = bar_chart_error(letters+['mnist'], list_classifiers, list_colors)
fig = go.Figure(data=[trace_bar[key] for key in trace_bar], layout=layout_bar)
plotly.offline.iplot(fig)
############# kNN
from scipy.spatial.distance import cdist
from sklearn.metrics import confusion_matrix
from collections import Counter
from matplotlib import cm
k_max = 5
x_train, y_train = data['mnist']['train']['X'], data['mnist']['train']['Y']
x_test, y_test = data['mnist']['test']['X'], data['mnist']['test']['Y']
def confusion_kNN(k, x_train=x_train, y_train=y_train,
x_test = x_test, y_test=y_test):
assert k>0
# fonction auxiliaire: prédiction du label (par vote majoritaire)
# en fonction des indices des plus proches voisins
def predicted_label(ind_neighbors):
# Conversion de l'array y_train[ind_neighbors] (contenant les
# labels des voisins d'entraînement) en tuple
label_neighbors = tuple(map(tuple, y_train[ind_neighbors]))
# élément majoritaire dans label_neighbors
return Counter(label_neighbors).most_common(1)[0][0]
# différences des distances 2 à 2 entre les images test et les images d'entraînement
dist_matrix = cdist(x_test, x_train, 'euclidean')
y_predicted = []
for i in range(len(x_test)):
# indices des k plus petites valeurs de dist_matrix[i, :]
# i.e. indices des k plus petites distances entre
# l'image-test d'indice i et les images d'entraînement
# i.e. indices des k images d'entraînement les plus proches
# de l'image-test i (par définition de dist_matrix)
ind_k_smallest = np.argpartition(dist_matrix[i, :], k)[:k]
y_predicted.append(predicted_label(ind_k_smallest))
return confusion_matrix(y_test, np.array(y_predicted))
# détermination de l'erreur de classification
def classification_err(k_max=k_max, x_train=x_train, y_train=y_train, x_test=x_test, y_test=y_test):
# liste des erreurs de classification
err = []
len_test = len(x_test)
for k in range(1, k_max+1):
# matrice de confusion à laquelle on enlève la diagonale
# (laquelle correspond aux prédictions correctes, dont la perte 0-1 vaut zéro)
confusion_err = confusion_kNN(k, x_train=x_train, y_train=y_train, x_test=x_test, y_test=y_test)
confusion_err -= np.diag(np.diag(confusion_err))
# calcul de l'erreur de prédiction
err.append(np.sum(confusion_err) / len_test)
return np.array(err)
errors_train = classification_err(x_test=x_train, y_test=y_train)
errors_test = classification_err()
# Diagramme en bâtons
def bar_chart_error_kNN(errors_train=errors_train, errors_test=errors_test, trace_bar = trace_bar):
learning_sets = ['mnist train', 'mnist test']
k_max = len(errors_train)
# Couleur des bâtons
cm = plt.get_cmap('Blues')
colors = [np.array(cm(1.*i/(k_max+1))) for i in range(2, k_max+3)]
rgb_colors = ['rgb{}'.format(tuple(255*c[:-1])) for c in colors]
for k, color, err in zip(range(k_max+1), rgb_colors, zip(errors_train, errors_test)):
trace_bar['{}NN'.format(k+1)] = go.Bar(
x = learning_sets,
y = list(err),
name = '{}NN'.format(k+1),
marker=dict(
color=color,
),
opacity=0.75
)
return trace_bar
trace_bar = bar_chart_error_kNN()
fig = go.Figure(data=[trace_bar[key] for key in trace_bar], layout=layout_bar)
plotly.offline.iplot(fig)
La méthode des $k$ plus proches voisins, qui est un modèle discriminatif puisqu'aucune supposition n'est faite sur la distribution des $x_i$, est plus performante sur cette tâche de classification.
En effet, les modèles génératifs QDA et LDA essayent d'estimer $P(X, Y)$ en faisant des suppostions les distributions qui générent (d'où le nom) les $x_i$, pour chaque classe (en l'occurrence, la supposition faite est qu'il s'agit de Gaussiennes).
En ne faisant pas de supposition sur la distribution des $x_i$, la méthode des $k$ plus proches voisins ne s'attache qu'au fait que plus des images sont proches (au sens de la distance euclidienne) plus elles sont susceptibles de représenter le même chiffre (même classe).
On notera que les modèles génératifs ne sont pas à leur avantage dans cette tâche de classification de chiffres, puisqu'il n'y a aucune raison que les images d'une même classe soient générées par une Gaussienne : par le théorème central-limite, la génération par des Gaussiennes est une hypothèse raisonnable quand les données résultent d'une multitude de processus aléatoires suivant la même loi (comme dans la nature par exemple), ce qui n'est pas le cas ici.