.. _tdnote20172rst: =============================== 1A.e - TD noté, 21 février 2017 =============================== .. only:: html **Links:** :download:`notebook `, :downloadlink:`html `, :download:`python `, :downloadlink:`slides `, :githublink:`GitHub|_doc/notebooks/exams/td_note_2017_2.ipynb|*` 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. .. code:: ipython3 from jyquickhelper import add_notebook_menu add_notebook_menu() .. contents:: :local: Précision : dans tout l’énoncé, la transposée d’une matrice est notée :math:`X' = X^{T}`. La plupart du temps :math:`X` et :math:`Y` désignent des vecteurs colonnes. :math:`\beta` désigne un vecteur ligne, :math:`W` une matrice diagonale. Exercice 1 ========== Q1 ~~ A l’aide du module `random `__, générer un ensemble de points aléatoires. .. code:: ipython3 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 .. parsed-literal:: [1000, 51, 83, 29, 15, 62, 90, 28, 61, 40] Q2 ~~ La médiane d’un ensemble de points :math:`\left\{X_1, ..., X_n\right\}` est une valeur :math:`X_M` telle que : .. math:: \sum_i \mathbf{1\!\!1}_{X_i < X_m} = \sum_i \mathbf{1\!\!1}_{X_i > X_m} Autrement dit, il y a autant de valeurs inférieures que supérieures à :math:`X_M`. On obtient cette valeur en triant les éléments par ordre croissant et en prenant celui du milieu. .. code:: ipython3 def mediane(ensemble): tri = list(sorted(ensemble)) return tri[len(tri)//2] mediane(ens) .. parsed-literal:: 61 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. .. code:: ipython3 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) .. parsed-literal:: 56.0 Q4 ~~ Pour un ensemble de points :math:`E=\left\{X_1, ..., X_n\right\}`, on considère la fonction suivante : .. math:: f(x) = \sum_{i=1}^n \left | x - X_i\right | . On suppose que la médiane :math:`X_M` de l’ensemble :math:`E` n’appartient pas à :math:`E` : :math:`X_M \notin E`. Que vaut :math:`f'(X_M)` ? On acceptera le fait que la médiane est le seul point dans ce cas. .. math:: f'(X_m) = - \sum_{i=1}^n \mathbf{1\!\!1}_{X_i < X_m} + \sum_{i=1}^n \mathbf{1\!\!1}_{X_i > X_m} Par définition de la médiane, :math:`f'(X_M)=0`. En triant les éléments, on montre que la :math:`f'(x) = 0 \Longleftrightarrow x=X_m`. Q5 ~~ On suppose qu’on dispose d’un ensemble d’observations :math:`\left(X_i, Y_i\right)` avec :math:`X_i, Y_i \in \mathbb{R}`. La régression linéaire consiste une relation linéaire :math:`Y_i = a X_i + b + \epsilon_i` qui minimise la variance du bruit. On pose : .. math:: E(a, b) = \sum_i \left(Y_i - (a X_i + b)\right)^2 On cherche :math:`a, b` tels que : .. math:: 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 : .. math:: \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 : .. math:: \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 : .. math:: \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 :math:`X`, on écrit le problème d’optimisation, on cherche les coefficients :math:`\beta^*` qui minimisent : .. math:: 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 : :math:`\beta^* = (X'X)^{-1}X'Y`. Ecrire une fonction qui calcule ce vecteur optimal. .. code:: ipython3 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 .. parsed-literal:: array([[ 1.00141843]]) Q6 ~~ Ecrire une fonction qui transforme un vecteur en une matrice diagonale. .. code:: ipython3 def matrice_diagonale(W): return numpy.diag(W) matrice_diagonale([1, 2, 3]) .. parsed-literal:: 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 :math:`w_i`. On veut maintenant trouver le vecteur :math:`\beta` qui minimise : .. math:: 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ù :math:`W=diag(w_1, ..., w_n)` est la matrice diagonale. La solution est : .. math:: \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. .. code:: ipython3 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) .. parsed-literal:: (array([[ 1.]]), array([[ 1.01240451]])) Q8 ~~ Ecrire une fonction qui calcule les quantités suivantes (fonctions `maximum `__, `reciprocal `__). .. math:: z_i = \frac{1}{\max\left( \delta, \left|y_i - X_i \beta\right|\right)} .. code:: ipython3 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) .. parsed-literal:: 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]]) Q9 ~~ On souhaite coder l’algorithme suivant : 1. :math:`w_i^{(1)} = 1` 2. :math:`\beta_{(t)} = (X'W^{(t)}X)^{-1}X'W^{(t)}Y` 3. :math:`w_i^{(t+1)} = \frac{1}{\max\left( \delta, \left|y_i - X_i \beta^{(t)}\right|\right)}` 4. :math:`t = t+1` 5. Retour à l’étape 2. .. code:: ipython3 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) .. parsed-literal:: 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]] .. parsed-literal:: array([[ 2.07611358]]) .. code:: ipython3 regression_lineaire(X, Y) .. parsed-literal:: array([[ 13.82243581]]) Q10 ~~~ .. code:: ipython3 ens = ensemble_aleatoire(10) Y = numpy.empty((len(ens), 1)) Y[:,0] = ens X = numpy.ones((len(ens), 1)) mediane(ens) .. parsed-literal:: 34.5 .. code:: ipython3 Y.mean(axis=0) .. parsed-literal:: array([ 131.1]) .. code:: ipython3 regression_lineaire(X, Y) .. parsed-literal:: array([[ 131.1]]) .. code:: ipython3 algorithm(X,Y) .. parsed-literal:: 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]] .. parsed-literal:: array([[ 39.7931558]]) .. code:: ipython3 mediane(ens) .. parsed-literal:: 34.5 .. code:: ipython3 list(sorted(ens)) .. parsed-literal:: [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. 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 :math:`X_M` qui minimise la quantité : .. math:: \sum_{i=1}^n \left| X_i - X_M \right | Nous reprenons l’algorithme décrit ci-dessus : 1. :math:`w_i^{(1)} = 1` 2. :math:`\beta_{(t)} = (X'W^{(t)}X)^{-1}X'W^{(t)}Y` 3. :math:`w_i^{(t+1)} = \frac{1}{\max\left( \delta, \left|y_i - X_i \beta^{(t)}\right|\right)}` 4. :math:`t = t+1` 5. Retour à l’étape 2. L’erreur quadratique pondéré est : .. math:: E_2(\beta, W) = \sum_{i=1}^n w_i \left\Vert Y_i - X_i \beta \right\Vert^2 Si :math:`w_i = \frac{1}{\left|y_i - X_i \beta\right|}`, on remarque que : .. math:: 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 :math:`\beta` qui minimise :math:`E_2(\beta, W^{(t)})`, par construction, il ressort que : .. math:: 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 :math:`t \rightarrow E_1(\beta^{(t)})` est suite décroissant est minorée par 0. Elle converge donc vers un minimum. Or la fonction :math:`\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 :math:`\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 [@](https://www.python.org/dev/peps/pep-0465/) a été introduit par Python 3.5 et est équivalent à la fonction `numpy.dot `__. Les dimensions des matrices posent souvent quelques problèmes. .. code:: ipython3 import numpy y = numpy.array([1, 2, 3]) M = numpy.array([[3, 4], [6, 7], [3, 3]]) M.shape, y.shape .. parsed-literal:: ((3, 2), (3,)) .. code:: ipython3 try: M @ y except Exception as e: print(e) .. parsed-literal:: 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 : .. code:: ipython3 y @ M .. parsed-literal:: array([24, 27]) Ou : .. code:: ipython3 M.T @ y .. parsed-literal:: array([24, 27])