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.
from jyquickhelper import add_notebook_menu
add_notebook_menu()
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.
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
[1000, 51, 83, 29, 15, 62, 90, 28, 61, 40]
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.
def mediane(ensemble):
tri = list(sorted(ensemble))
return tri[len(tri)//2]
mediane(ens)
61
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.
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)
56.0
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$.
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.
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
array([[ 1.00141843]])
Ecrire une fonction qui transforme un vecteur en une matrice diagonale.
def matrice_diagonale(W):
return numpy.diag(W)
matrice_diagonale([1, 2, 3])
array([[1, 0, 0], [0, 2, 0], [0, 0, 3]])
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.
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)
(array([[ 1.]]), array([[ 1.01240451]]))
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)}$$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)
array([[ 1.01330469e-03], [ 5.26315789e+00], [ 4.54545455e+00], [ 3.22580645e+00], [ 1.85185185e+00], [ 1.47058824e+00], [ 1.14942529e+00], [ 1.07526882e+00], [ 1.07526882e+00], [ 1.00000000e-01]])
On souhaite coder l'algorithme suivant :
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 150.13052808 [[ 13.82243581]] 1 104.79608014 [[ 3.21524459]] 2 100.851019446 [[ 2.25815451]] 3 100.36420567 [[ 2.12644545]] 4 100.255554539 [[ 2.09141327]] 5 100.220626093 [[ 2.0777948]] 6 100.219023635 [[ 2.07639404]] 7 100.21901041 [[ 2.07631459]] 8 100.218994922 [[ 2.07622156]] 9 100.218976948 [[ 2.07611358]]
array([[ 2.07611358]])
regression_lineaire(X, Y)
array([[ 13.82243581]])
ens = ensemble_aleatoire(10)
Y = numpy.empty((len(ens), 1))
Y[:,0] = ens
X = numpy.ones((len(ens), 1))
mediane(ens)
34.5
Y.mean(axis=0)
array([ 131.1])
regression_lineaire(X, Y)
array([[ 131.1]])
algorithm(X,Y)
0 1737.8 [[ 131.1]] 1 1215.2110733 [[ 55.05276833]] 2 1196.55478823 [[ 48.77739411]] 3 1190.4919578 [[ 45.7459789]] 4 1183.56462833 [[ 42.28231416]] 5 1179.0 [[ 39.7931558]] 6 1179.0 [[ 39.7931558]] 7 1179.0 [[ 39.7931558]] 8 1179.0 [[ 39.7931558]] 9 1179.0 [[ 39.7931558]]
array([[ 39.7931558]])
mediane(ens)
34.5
list(sorted(ens))
[5, 6, 12, 14, 29, 40, 52, 67, 86, 1000]
La régression linéaire égale la moyenne, l'algorithme s'approche de la médiane.
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 :
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.
import numpy
y = numpy.array([1, 2, 3])
M = numpy.array([[3, 4], [6, 7], [3, 3]])
M.shape, y.shape
((3, 2), (3,))
try:
M @ y
except Exception as e:
print(e)
shapes (3,2) and (3,) not aligned: 2 (dim 1) != 3 (dim 0)
Par défaut, numpy considère un vecteur de taille (3,)
comme un vecteur ligne (3,1)
. Donc l'expression suivante va marcher :
y @ M
array([24, 27])
Ou :
M.T @ y
array([24, 27])