1A.e - TD noté, 21 février 2017¶

Solution du TD noté, celui-ci présente un algorithme pour calculer les coefficients d’une régression quantile et par extension d’une médiane dans un espace à plusieurs dimensions.

Précision : dans tout l’énoncé, la transposée d’une matrice est notée $X' = X^{T}$ . La plupart du temps $X$ et $Y$ désignent des vecteurs colonnes. $\beta$ désigne un vecteur ligne, $W$ une matrice diagonale.

Exercice 1¶

Q1¶

A l’aide du module random, générer un ensemble de points aléatoires.

[1]:

import random


def ensemble_aleatoire(n):
    res = [random.randint(0, 100) for i in range(n)]
    res[0] = 1000
    return res


ens = ensemble_aleatoire(10)
ens

[1]:

[1000, 15, 65, 44, 95, 17, 43, 44, 88, 88]

Q2¶

La médiane d’un ensemble de points $\left\{X_1, ..., X_n\right\}$ est une valeur $X_M$ telle que :

$\sum_i \mathbb{1}_{X_i < X_m} = \sum_i \mathbb{1}_{X_i > X_m}$

Autrement dit, il y a autant de valeurs inférieures que supérieures à $X_M$ . On obtient cette valeur en triant les éléments par ordre croissant et en prenant celui du milieu.

[2]:

def mediane(ensemble):
    tri = list(sorted(ensemble))
    return tri[len(tri) // 2]


mediane(ens)

[2]:

Q3¶

Lorsque le nombre de points est pair, la médiane peut être n’importe quelle valeur dans un intervalle. Modifier votre fonction de façon à ce que la fonction précédente retourne le milieu de la fonction.

[3]:

def mediane(ensemble):
    tri = list(sorted(ensemble))
    if len(tri) % 2 == 0:
        m = len(tri) // 2
        return (tri[m] + tri[m - 1]) / 2
    else:
        return tri[len(tri) // 2]


mediane(ens)

[3]:

54.5

Q4¶

Pour un ensemble de points $E=\left\{X_1, ..., X_n\right\}$ , on considère la fonction suivante :

$f(x) = \sum_{i=1}^n \left | x - X_i\right |$

On suppose que la médiane $X_M$ de l’ensemble $E$ n’appartient pas à $E$ : $X_M \notin E$ . Que vaut $f'(X_M)$ ? On acceptera le fait que la médiane est le seul point dans ce cas.

$f'(X_m) = - \sum_{i=1}^n \mathbb{1}_{X_i < X_m} + \sum_{i=1}^n \mathbb{1}_{X_i > X_m}$

Par définition de la médiane, $f'(X_M)=0$ . En triant les éléments, on montre que la $f'(x) = 0 \Longleftrightarrow x=X_m$ .

Q5¶

On suppose qu’on dispose d’un ensemble d’observations $\left(X_i, Y_i\right)$ 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 \left(Y_i - (a X_i + b)\right)^2$

On cherche $a, b$ tels que :

$a^*, b^* = \arg \min E(a, b) = \arg \min \sum_i \left(Y_i - (a X_i + b)\right)^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}$

Lorsqu’on a plusieurs dimensions pour $X$ , on écrit le problème d’optimisation, on cherche les coefficients $\beta^*$ qui minimisent :

$E(\beta)=\sum_{i=1}^n \left(y_i - X_i \beta\right)^2 = \left \Vert Y - X\beta \right \Vert ^2$

La solution est : $\beta^* = (X'X)^{-1}X'Y$ .

Ecrire une fonction qui calcule ce vecteur optimal.

[4]:

from numpy.linalg import inv


def regression_lineaire(X, Y):
    t = X.T
    return inv(t @ X) @ t @ Y


import numpy

X = numpy.array(ens).reshape((len(ens), 1))
regression_lineaire(
    X, X + 1
)  # un essai pour vérifier que la valeur n'est pas aberrante

[4]:

array([[1.00144835]])

Q6¶

Ecrire une fonction qui transforme un vecteur en une matrice diagonale.

[5]:

def matrice_diagonale(W):
    return numpy.diag(W)


matrice_diagonale([1, 2, 3])

[5]:

array([[1, 0, 0],
       [0, 2, 0],
       [0, 0, 3]])

Q7¶

On considère maintenant que chaque observation est pondérée par un poids $w_i$ . On veut maintenant trouver le vecteur $\beta$ qui minimise :

$E(\beta)=\sum_{i=1}^n w_i \left( y_i - X_i \beta \right)^2 = \left \Vert W^{\frac{1}{2}}(Y - X\beta)\right \Vert^2$

Où $W=diag(w_1, ..., w_n)$ est la matrice diagonale. La solution est :

$\beta_* = (X'WX)^{-1}X'WY$

Ecrire une fonction qui calcule la solution de la régression pondérée. La fonction ravel est utile.

[6]:

def regression_lineaire_ponderee(X, Y, W):
    if len(W.shape) == 1 or W.shape[0] != W.shape[1]:
        # c'est un vecteur
        W = matrice_diagonale(W.ravel())
    wx = W @ X
    xt = X.T
    return inv(xt @ wx) @ xt @ W @ Y


X = numpy.array(sorted(ens)).reshape((len(ens), 1))
Y = X.copy()
Y[0] = max(X)
W = numpy.ones(len(ens))
W[0] = 0
regression_lineaire_ponderee(X, Y, W), regression_lineaire(X, Y)

[6]:

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

Q8¶

Ecrire une fonction qui calcule les quantités suivantes (fonctions maximum, reciprocal).

$z_i = \frac{1}{\max\left( \delta, \left|y_i - X_i \beta\right|\right)}$

[7]:

def calcule_z(X, beta, Y, W, delta=0.0001):
    epsilon = numpy.abs(Y - X @ beta)
    return numpy.reciprocal(numpy.maximum(epsilon, numpy.ones(epsilon.shape) * delta))


calcule_z(X * 1.0, numpy.array([[1.01]]), Y, W)

[7]:

array([[1.01538305e-03],
       [5.88235294e+00],
       [2.32558140e+00],
       [2.27272727e+00],
       [2.27272727e+00],
       [1.53846154e+00],
       [1.13636364e+00],
       [1.13636364e+00],
       [1.05263158e+00],
       [1.00000000e-01]])

Q9¶

On souhaite coder l’algorithme suivant :

$w_i^{(1)} = 1$
$\beta_{(t)} = (X'W^{(t)}X)^{-1}X'W^{(t)}Y$
$w_i^{(t+1)} = \frac{1}{\max\left( \delta, \left|y_i - X_i \beta^{(t)}\right|\right)}$
$t = t+1$
Retour à l’étape 2.

[8]:

def algorithm(X, Y, delta=0.0001):
    W = numpy.ones(X.shape[0])
    for i in range(0, 10):
        beta = regression_lineaire_ponderee(X, Y, W)
        W = calcule_z(X, beta, Y, W, delta=delta)
        E = numpy.abs(Y - X @ beta).sum()
        print(i, E, beta)
    return beta


X = numpy.random.rand(10, 1)
Y = X * 2 + numpy.random.rand()
Y[0] = Y[0] + 100
algorithm(X, Y)

0 173.8971776057001 [[24.09487215]]
1 113.8229016016727 [[6.16379991]]
2 102.31983967417062 [[2.73034639]]
3 100.39132855434099 [[2.15204012]]
4 100.24844521352179 [[2.10443138]]
5 100.10481587527799 [[2.05258271]]
6 100.0856692693878 [[2.04297748]]
7 100.08332794594425 [[2.04131809]]
8 100.0800296132712 [[2.03867545]]
9 100.07493795090468 [[2.03459599]]

[8]:

array([[2.03459599]])

[9]:

regression_lineaire(X, Y)

[9]:

array([[24.09487215]])

Q10¶

[10]:

ens = ensemble_aleatoire(10)
Y = numpy.empty((len(ens), 1))
Y[:, 0] = ens
X = numpy.ones((len(ens), 1))
mediane(ens)

[10]:

51.0

[11]:

Y.mean(axis=0)

[11]:

array([136.7])

[12]:

regression_lineaire(X, Y)

[12]:

array([[136.7]])

[13]:

algorithm(X, Y)

0 1726.6000000000004 [[136.7]]
1 1168.2374836735025 [[58.61874184]]
2 1163.7683647305184 [[56.38418237]]
3 1161.0 [[54.75688683]]
4 1161.0 [[54.75688683]]
5 1161.0 [[54.75688683]]
6 1161.0 [[54.75688683]]
7 1161.0 [[54.75688683]]
8 1161.0 [[54.75688683]]
9 1161.0 [[54.75688683]]

[13]:

array([[54.75688683]])

[14]:

mediane(ens)

[14]:

51.0

[15]:

list(sorted(ens))

[15]:

[3, 7, 22, 24, 47, 55, 62, 70, 77, 1000]

La régression linéaire égale la moyenne, l’algorithme s’approche de la médiane.

Quelques explications et démonstrations¶

Cet énoncé est inspiré de Iteratively reweighted least squares. Cet algorithme permet notamment d’étendre la notion de médiane à des espaces vectoriels de plusieurs dimensions. On peut détermine un point $X_M$ qui minimise la quantité :

$\sum_{i=1}^n \left| X_i - X_M \right |$

Nous reprenons l’algorithme décrit ci-dessus :

$w_i^{(1)} = 1$
$\beta_{(t)} = (X'W^{(t)}X)^{-1}X'W^{(t)}Y$
$w_i^{(t+1)} = \frac{1}{\max\left( \delta, \left|y_i - X_i \beta^{(t)}\right|\right)}$
$t = t+1$
Retour à l’étape 2.

L’erreur quadratique pondéré est :

$E_2(\beta, W) = \sum_{i=1}^n w_i \left\Vert Y_i - X_i \beta \right\Vert^2$

Si $w_i = \frac{1}{\left|y_i - X_i \beta\right|}$ , on remarque que :

$E_2(\beta, W) = \sum_{i=1}^n \frac{\left\Vert Y_i - X_i \beta \right\Vert^2}{\left|y_i - X_i \beta\right|} = \sum_{i=1}^n \left|y_i - X_i \beta\right| = E_1(\beta)$

On retombe sur l’erreur en valeur absolue optimisée par la régression quantile. Comme l’étape 2 consiste à trouver les coefficients $\beta$ qui minimise $E_2(\beta, W^{(t)})$ , par construction, il ressort que :

$E_1(\beta^{(t+1)}) = E_2(\beta^{(t+1)}, W^{(t)}) \leqslant E_2(\beta^{(t)}, W^{(t)}) = E_1(\beta^{(t)})$

La suite $t \rightarrow E_1(\beta^{(t)})$ est suite décroissant est minorée par 0. Elle converge donc vers un minimum. Or la fonction $\beta \rightarrow E_1(\beta)$ est une fonction convexe. Elle n’admet qu’un seul minimum (mais pas nécessaire un seul point atteignant ce minimum). L’algorithme converge donc vers la médiane. Le paramètre $\delta$ est là pour éviter les erreurs de divisions par zéros et les approximations de calcul faites par l’ordinateur.

Quelques commentaires sur le code¶

Le symbol @ a été introduit par Python 3.5 et est équivalent à la fonction numpy.dot. Les dimensions des matrices posent souvent quelques problèmes.

[16]:

import numpy

y = numpy.array([1, 2, 3])
M = numpy.array([[3, 4], [6, 7], [3, 3]])
M.shape, y.shape

[16]:

((3, 2), (3,))

[17]:

try:
    M @ y
except Exception as e:
    print(e)

matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 3 is different from 2)

Par défaut, numpy considère un vecteur de taille (3,) comme un vecteur ligne (3,1). Donc l’expression suivante va marcher :

[18]:

y @ M

[18]:

array([24, 27])

Ou :

[19]:

M.T @ y

[19]:

array([24, 27])

[ ]:

Notebook on github