La machine à vecteurs de support

add.png

\bullet Présentation:

La Machine à Vecteurs de Supports (SVM), également connue sous le nom de séparateurs à vaste marge, a été élaborée par Vladimir Vapnik en 1995. Elle est utilisée pour la classification d’une variable réponse binaire Y à partir d’une matrice \mathbf{X} = (X ^1, \cdots, X ^P) de P variables explicatives continues.

Les SVM s’appuient sur la notion de marge maximale et la notion de fonction noyau (kernel), permettant de traiter des problèmes de discrimination non-linéaire en reformulant le problème de classement comme un problème d’optimisation quadratique.

Les SVM font parties des techniques d’apprentissage supervisé et nécessitent le passage par une méthode de validation afin de généraliser leur résultat. En effet, il s’agit d’un outil très puissant permettant de modéliser avec une extrême précision un phénomène en fonction du noyau considéré, augmentant ainsi fortement le risque de surapprentissage.

\bullet Les machines à vecteurs de support:

Hypothèse préliminaire: Y binaire, \mathbf{X} continue.

L’Algorithme:

Le calcul des coefficients via les SVM est particulièrement coûteux étant donné qu’il se base sur la recherche d’une solution optimale au travers de noyau pouvant rallonger les ressources de calcul demandée, notamment pour le noyau polynomial. De plus, l’idée principale des SVM est d’estimer les vecteurs supports qui représentent les pondérations associées à chaque observation afin de déterminer leur influence sur l’élaboration de la frontière décisionnelle. En outre, plus il y a d’observations à considérer et plus le système associé à l’estimation de ces pondérations est important.

Quatre étapes sont à suivre pour estimer les coefficients \alpha_i ^*, i \in [1,n^*], qui ne correspondent donc pas aux P variables mais à nos n^* < n individus pour lesquels l’influence est significative (ceux associés à un vecteur support négligeable sont alors supprimés du modèle prédictif final).

– Étape 0: Définir le noyau et les paramètres qui lui sont associés.

– Étape 1: Déterminer l’équation des vecteurs de support à partir des paramètres définis en étape 0,

Q_{\alpha} = \sum_{i = 1} ^n \alpha_i - \frac{1}{2} \sum_{i_1 = 1} ^n \sum_{i_2 = 1} ^n \alpha_{i_1} \alpha_{i_2} Y_{i_1} Y_{i_2} K(X_{i_1}, X_{i_2})

– Étape 2: Déterminer la matrice à résoudre et dont la ième ligne est de la forme,

\alpha_i ^* = \frac{Q_{\alpha}}{\partial Q_{\alpha_i}} = \alpha_i - \frac{1}{2} \sum_{i_2 = 1} ^n \alpha_{i_2} Y_i Y_{i_2} K(X_i, X_{i_2}) = 1, \forall i \in [1,n]

, où \alpha = (\alpha_1, \cdots, \alpha_n) représente le vecteur des multiplicateurs de Lagrange.

– Étape 3: Rajouter à la matrice, conçue en étape 2, le vecteur en ligne et basé sur la contrainte de construction des \alpha, \sum_{i = 1} ^n Y_i \alpha_i = 0. Ce système est alors noté \blacksquare.

– Étape 4: Résoudre le système \blacksquare via les multiplicateurs de Lagrange afin d’obtenir les estimations (\alpha_1 ^*, \cdots, \alpha_n ^*).

Équation prédictive:

L’Algorithme permet de déterminer les Lagrangiens \alpha ^* soit les coefficients d’influence de nos vecteurs supports associés aux différentes observations. Cependant, les estimations \alpha ^* ne sont pas suffisants pour la prédiction d’une nouvelle observation. Définissons E ^* l’ensemble des observations associées à un vecteur support non nul et dont l’influence doit être considérée pour la prédiction. L’Équation de prédiction est alors,

f(x) = \sum_{i \in E ^*} \alpha_i ^* Y_i < X_i, x> + \beta_0

Où,

x est la notation utilisée pour désigner les valeurs caractéristiques associées aux variables explicatives considérées pour la nouvelle observation à prédire,

<.,.> le produit scalaire,

\beta_0 le coefficient constant obtenu par,

i; \beta_0 = Y_i - \sum_{p = 1} ^P X_i ^p \beta_p

, et dont la valeur ne varie pas quelque soit l’observation i utilisée.

\beta = (\beta_1, \cdots, \beta_P) les coefficients obtenus par,

\forall p \in [1,P], \beta_p = \sum_{i = 1} ^n \alpha_i ^* X_i ^p

Le choix du noyau:

Deux cas sont à distinguer lors de l’usage des SVM: le cas linéairement séparable et non linéairement séparable. Pour le premier, cela revient à construire une frontière linéaire optimale selon la notion de marge maximale. Le noyau est alors le noyau linéaire de formulation,

K(X_{i_1},X_{i_2}) = X{i_1} \cdot X_{i_2}

Pour le second cas, l’idée est de passer la construction d’une nouvelle dimension qui va envoyer les données dans un espace où la séparation est possible. Trois noyaux principaux existent parmi la longue liste disponible,

– le noyau polynomial K(X_{i_1},X_{i_2}) = (\gamma \cdot X_{i_1} ^T \cdot X_{i_2} + w_0) ^d,

– le noyau radial K(X_{i_1},X_{i_2}) = e^{- \gamma \cdot \Vert X_{i_1} - X_{i_2} \Vert ^2}, avec,

\Vert X_{i_1} - X_{i_2} \Vert = \sqrt{<X_{i_1},X_{i_1}> + <X_{i_2},X_{i_2}> - 2 <X_{i_1},X_{i_2}>}

– le noyau sigmoïdale K(X_{i_1},X_{i_2}) = tanh(\gamma \cdot X_{i_1} ^T \cdot X_{i_2} + w_0).

\gamma, w_0 et d sont les paramètres à figer.

La sélection de variable:

Le choix des variables à retenir dans le modèle final amène à l’introduction du principe de « marge ». La marge correspond à la distance entre la frontière construite par le modèle et ses vecteurs supports. En général, un modèle est de bonne qualité si la marge est suffisamment éloignée des vecteurs supports, permettant ainsi plus d’aisance dans sa généralisation.

La sélection des variables se fait en fonction de la marge associée et notamment de sa sensibilité. L’Idée est de s’orienter vers un modèle dont la marge est la plus robuste possible.

La marge se calcul selon la formule suivante,

m = \frac{2}{|| \beta ||} = \frac{2}{\sqrt{\sum_{i_1 = 1} ^n \sum_{i_2 = 1} ^n \alpha_{i_1} \alpha_{i_2} Y_{i_1} Y_{i_2} <X_{i_1},X_{i_2}>}}

En définissant un facteur échelle v qui est en fait le vecteur unitaire pondérant le noyau utilisé K, le système à résoudre pour définir la sensibilité de la marge est alors,

| \frac{\partial || w || ^2}{\partial v} | = | -2 \sum_{i_1, i_2} \alpha_{i_1} ^* \alpha_{i_2} ^* Y_{i_1} Y_{i_2} \frac{\partial K(v X_{i_1},v X_{i_2})}{\partial v} |

Généralisation au cas multiclasse:

Les SVM ne sont pas directement transposables au cas multiclasse du fait d’une formule d’estimation basée sur un principe binaire. Pour se faire, l’idée est d’appliquer une approche du type « 1-contre-1 » et qui consiste, pour K classes, à construire \frac{K \times (K - 1)}{2} classifieurs et d’avoir recours à un schéma de vote afin de décider de la classe de prédiction.

\bullet Annexe théorique:

Nous présentons ici une esquisse de la démarche méthodologique conduisant à l’algorithme d’estimation des coefficients \alpha d’un modèle basé sur les SVM.

Nous nous plaçons dans l’espace E munit du produit scalaire. L’objectif est de construire une fonction h qui sépare au mieux \mathbf{X}.

Le cas linéairement séparable:

Le plus simple, il s’écrit sous la forme suivante,

h(\mathbf{X}) = w ^t \cdot \mathbf{X} + w_0

Avec,

h(X_i) \geq 0 \Rightarrow Y_i = 1,

h(X_i) < 0 \Rightarrow Y_i = -1,

La frontière de décision h(\mathbf{X}) = 0 est un hyperplan, appelé hyperplan séparateur ou séparatrice. L’objectif est d’avoir,

\forall i \in [1,n], y_i \cdot h(X_i) \geq 0

Notion de marge maximale:

L’équation du séparateur linéaire montre que nous pouvons générer une infinité de classifieurs, la notion de marge prend alors toute son importance. Elle se définit comme étant la distance entre l’hyperplan et les observations les plus proches qui sont en fait nos vecteurs supports. L’hyperplan, à déterminer et qui maximise la marge, est donné par,

arg \hspace*{1mm} max_{w,w_0} min_i {\parallel \mathbf{X} - X_i \parallel : \mathbf{X} \in \mathbb{R} ^n / h(\mathbf{X}) = w^t \cdot \mathbf{X} + w_0 = 0}

Formulation primale:

La distance de X_i à l’hyperplan est donnée par sa projection orthogonale sur l’hyperplan,

\frac{y_i \cdot (w^T \cdot X_i + w_0)}{\parallel w \parallel}

La formule que nous recherchons se construit sous la contrainte:

arg \hspace*{1mm} max_{w,w_0} {\frac{1}{\parallel w \parallel} min_i [y_i \cdot (w^t \cdot X_i + w_0)]}

Afin de faciliter l’optimisation, nous choisissons de normaliser w et w_0, de telle manière à ce que les observations à la marge (\mathbf{X}_+ pour les vecteurs supports sur la frontière positive, et \mathbf{X}_- pour ceux situés sur la frontière opposée) satisfassent,

w^t \cdot \mathbf{X}_+ + w_0 = +1

w^t \cdot \mathbf{X}_- + w_0 = -1

La marge à maximiser vaut alors \frac{1}{\parallel w \parallel}. La formulation dite primale des SVM est ainsi,

Minimiser \frac{1}{2} \parallel w \parallel ^2 sous contraintes que y_i \cdot (w^t \cdot X_i + w_0) \geq 1

Ceci peut se résoudre par la méthode classique des multiplicateurs de Lagrange \alpha, où le lagrangien est donné par,

L(w,w_0,\alpha) = \frac{1}{2} \parallel w \parallel ^2 - \sum_{i = 1} ^n \alpha_i \cdot {y_i \cdot (w^t  \cdot X_i + w_0) - 1} \hspace*{10mm} (\blacksquare)

Formulation duale:

Le lagrangien doit être minimisé par rapport à w et w_0, et maximisé par rapport à \alpha. Pour cela nous posons notre problème sous forme duale. Nous cherchons ainsi à annuler les dérivées partielles du lagrangien, selon les conditions de Kuhn-Tucker, Nous obtenons,

\sum_{i = 1} ^n \alpha_i \cdot y_i \cdot X_i = w^*

\sum_{i = 1} ^n \alpha_i \cdot y_i = 0

En réinjectant ces valeurs dans l’équation (\blacksquare), nous avons la formulation duale suivante,

Maximiser L(\alpha) = \sum_{i = 1} ^n \alpha_i - \frac{1}{2} \sum_{i_1,i_2} ^n \alpha_{i_1} \cdot \alpha_{i_2} \cdot y_{i_1} \cdot y_{i_2} \cdot X_{i_1} ^t \cdot X_{i_2}

, sous les contraintes \alpha_i \geq 0, et \sum_{i = 1} ^n \alpha_i y_i.

Ce qui nous donne les multiplicateurs de Lagrange optimaux \alpha_i ^*. L’Equation de l’hyperplan solution devient alors,

h(x) = \sum_{i = 1} ^n \alpha_i ^* \cdot y_i \cdot (x \cdot X_i) + w_0

Le cas non linéairement séparable:

Finalement l’hyperplan solution ne dépend que du produit scalaire entre \mathbf{X} et les vecteurs supports. Cette remarque est à l’origine de la deuxième innovation majeure des SVM: l’utilisation de la fonction noyau pour projeter \mathbf{X} dans un espace dit de redescription.

Plus formellement, nous appliquons aux vecteurs d’entrée \mathbf{X}, une transformation non-linéaire notée \varphi.

Dans cet espace, nous cherchons alors l’hyperplan,

h(\mathbf{X}) = w^t \cdot \varphi(\mathbf{X}) + w_0

, qui vérifie \forall i \in [1,n], y_i \cdot h(X_i) > 0.

En utilisant la même procédure que dans le cas sans transformation, nous aboutissons au problème d’optimisation suivant,

Maximiser L(\alpha) = \sum_{i = 1} ^n \alpha_i - \frac{1}{2} \sum_{i_1,i_2} ^n \alpha_{i_1} \cdot \alpha_{i_2} \cdot y_{i_1} \cdot y_{i_2} \cdot \varphi(X_{i_1})^t \cdot \varphi(X_{i_2})

, sous contraintes \alpha_i \geq 0, et \sum_{i = 1} ^n \alpha_i \cdot y_i = 0.

La fonction noyau:

Le problème de la formulation ci-dessus est qu’elle implique un produit scalaire entre vecteurs dans l’espace de redescription et de dimension élevée, ce qui est couteux en termes de calculs. Pour résoudre ce problème, nous utilisons une astuce connue sous le nom de « Kernel trick » consistant à utiliser une fonction noyau vérifiant,

K(X_{i_1}, X_{i_2}) = \varphi(X_{i_1})^t \cdot \varphi(X_{i_2})

D’où l’expression de l’hyperplan séparateur en fonction de la fonction noyau,

h(x) = \sum_{i = 1} ^n \alpha_i ^* \cdot y_i \cdot K(X_i, x) + w_0

L’intérêt de la fonction noyau est double,

– le calcul se fait dans l’espace d’origine, ceci est beaucoup moins coûteux qu’un produit scalaire en grande dimension,

– la transformation \varphi n’a pas besoin d’être connue explicitement, seule la fonction noyau intervient dans les calculs, augmentant nettement les possibilités de discrimination.

Enfin, en accords avec le théorème de Mercer, \varphi doit être symétrique et semi-définie positive. L’approche par Kernel trick généralise ainsi l’approche linéaire.

Extension aux marges souples:

Cependant, il n’est pas non plus toujours possible de trouver une séparatrice linéaire dans l’espace de redescription. Il se peut aussi que des observations soient mal étiquetées et que l’hyperplan séparateur ne soit pas la meilleure solution au problème de classement.

Nous nous ramenons alors à la technique dite de marge souple, qui tolère les mauvais classements. La méthode cherche un hyperplan séparateur qui minimise le nombre d’erreurs grâce à l’introduction de variables ressorts \xi_i, qui permettent de relâcher les contraintes sur les vecteurs d’apprentissage:

\forall \xi_i \geq 0, i \in [1,n], y_i \cdot (w^t \cdot X_i + w_0) \geq 1 - \xi_i

Avec les contraintes précédentes, le problème d’optimisation est modifié par un terme de pénalité qui agit directement sur les variables ressorts élevées et de formulation:

Minimiser \frac{1}{2} \parallel w \parallel ^2 + C \cdot \sum_{i = 1} ^n \xi_i, C > 0

, où C est une constante qui permet de contrôler le compromis entre nombre d’erreurs de classement et la largeur de la marge.

\bullet Exemple:

Soit le jeu de données suivant,

add

Certains le reconnaîtront, il s’agit de la fonction XOR également appelée fonction OU exclusif. Essayons de résoudre le problème en fonction des trois fonctions noyaux décrites pour le cas non linéairement séparable. La matrice de données qui lui est associée est la suivante,

\mathbf{X} = \begin{pmatrix} Y & X1 & X2 \\ -1 & -1 & -1 \\ 1 & -1 & 1 \\ 1 & 1 & -1 \\ -1 & 1 & 1 \\ \end{pmatrix}

Cas du noyau sigmoïdale:

Soit le noyau sigmoïdale de constante C = 1 et de paramètre \gamma = 1,

K(X_{i_1}, X_{i_2}) = tanh(1 + X_{i_1} ^1 X_{i_2} ^1 + X_{i_1} ^2 X_{i_2} ^2)

Nous calculons pour chaque combinaison d’observations le valeur du noyau,

\mathbf{Q} = \begin{pmatrix} i_1 & i_2 & Y_{i_1} \times Y_{i_2} \times K(X_{i_1}, X_{i_2}) \\ 1 & 1 & 0.9950548 \\ 1 & 2 & -0.7615942 \\ 1 & 3 & -0.76159421 \\ 1 & 4 & -0.7615942 \\ 2 & 1 & -0.7615942 \\ 2 & 2 & 0.9950548 \\ 2 & 3 & -0.7615942 \\ 2 & 4 & -0.7615942 \\ 3 & 1 & -0.7615942 \\ 3 & 2 & -0.7615942 \\ 3 & 3 & 0.9950548 \\ 3 & 4 & -0.7615942 \\ 4 & 1 & -0.7615942 \\ 4 & 2 & -0.7615942 \\ 4 & 3 & -0.7615942 \\ 4 & 4 & 0.9950548 \\ \end{pmatrix}

En rappelant la formule générale des Q_{\alpha},

Q_{\alpha} = \sum_{i = 1} ^n \alpha_i - \frac{1}{2} \cdot \sum_{i_1, i_2} ^n \alpha_{i_1} \alpha_{i_2} Q(i_1, i_2)

, nous obtenons la dérivation suivante,

\alpha_i ^* = \frac{\partial Q_{\alpha}}{\partial \alpha_i} = 1

\Rightarrow \alpha_i ^* = 1 - \sum_{i_2 = 1} ^n \alpha_{i_2} \times Q(i,i_2) = 1

\Rightarrow \alpha_i ^* = - \sum_{i_2 = 1} ^n \alpha_{i_2} \times Q(i,i_2)

Si nous raisonnons en termes de matrice afin de simplifier les calculs, cela revient à transposer les éléments de \mathbf{Q} en fonction de i_1. Nous obtenons alors la matrice,

\begin{pmatrix} \alpha_1 & \alpha_2 & \alpha_3 & \alpha_4 \\ -0.9950548 & 0.7615942 & 0.7615942 & 0.7615942 \\ 0.7615942 & -0.9950548 & 0.7615942 & 0.7615942 \\ 0.7615942 & 0.7615942 & -0.9950548 & 0.7615942 \\ 0.7615942 & 0.7615942 & 0.7615942 & -0.9950548 \\ \end{pmatrix}

Comme sous sommes sous contraintes que \sum_{i = 1} ^5 Y_i \times \alpha_i ^* = 0, nous rajoutons la condition au système à résoudre ainsi que la valeur de l’égalité (qui vaut -1 car, pour respecter le passage des équations à la forme matricielle du problème, nous faisons pivoter le coefficient constant associé à la somme des \alpha_i à droite),

\mathbf{S} = \begin{pmatrix} Y_t & \alpha_1 & \alpha_2 & \alpha_3 & \alpha_4 \\ -1 & -0.9950548 & 0.7615942 & 0.7615942 & 0.7615942 \\ -1 & 0.7615942 & -0.9950548 & 0.7615942 & 0.7615942 \\ -1 & 0.7615942 & 0.7615942 & -0.9950548 & 0.7615942 \\ -1 & 0.7615942 & 0.7615942 & 0.7615942 & -0.9950548 \\ 0 & 1 & -1 & -1 & 1 \end{pmatrix}

Nous résolvons ce système avec une régression linéaire et obtenons,

\alpha_1 ^* = \alpha_2 ^* = \alpha_3 ^* = \alpha_4 ^* = -0.7754

Tous nos coefficients sont différents de 0, nous en déduisons que pour ce noyau les quatre observations jouent un rôle de vecteur support.

Maintenant que les coefficients ont été estimés, nous pouvons déterminer la valeur de prédiction de X_1. Nous calculons dans un premier temps les coefficients \beta associés aux \alpha ^*,

\beta = (-0.7754 \times (-1) \times (-1) + (-0.7754) \times 1 \times (-1) + (-0.7754) \times 1 \times 1 + (-0.7754) \times (-1) \times 1, -0.7754 \times (-1) \times (-1) + (-0.7754) \times 1 \times 1 + (-0.7754) \times 1 \times (-1) + (-0.7754) \times (-1) \times 1)

= (0,0)

\Rightarrow \beta_0 = Y_1 - X_1 ^t \times \beta = Y_1 = -1

Nous pouvons maintenant déterminer la prédiction de X ^1 sur la base du modèle construit,

f(X_1) = \sum_{i = 1} ^4 \alpha_i ^* Y_i <X_i,X_1> + \beta_0

= -0.7754 \times (-1) \times ((-1) ^2 + (-1) ^2) + (-0.7754) \times 1 \times (-1 \times (-1) + (-1) \times 1) + (-0.7754) \times 1 \times ((-1) \times 1 + (-1) \times (-1)) + (-0.7754) \times (-1) \times ((-1) \times 1 + (-1) \times 1) + (-1)

= 0.7754 \times 2 - 0.7754 \times 0 -0.7754 \times 0 + 0.7754 \times (-2) - 1

= 1.5508 - 1.5508 - 1

\Rightarrow f(X_1) = - 1

Cas du noyau radial:

Soit le noyau radial de paramètre \gamma = 1,

K(X_{i_1}, X_{i_2}) = e ^{|| X_{i_1} - X_{i_2} || ^2}

Nous calculons pour chaque combinaison d’observations le valeur du noyau,

\mathbf{Q} = \begin{pmatrix} i_1 & i_2 & Y_{i_1} \times Y_{i_2} \times K(X_{i_1}, X_{i_2}) \\ 1 & 1 & 1 \\ 1 & 2 & -0.0183156389 \\ 1 & 3 & -0.0183156389 \\ 1 & 4 & 0.0003354626 \\ 2 & 1 & -0.0183156389 \\ 2 & 2 & 1 \\ 2 & 3 & 0.0003354626 \\ 2 & 4 & -0.0183156389 \\ 3 & 1 & -0.0183156389 \\ 3 & 2 & 0.0003354626 \\ 3 & 3 & 1 \\ 3 & 4 & -0.0183156389 \\ 4 & 1 & 0.0003354626 \\ 4 & 2 & -0.0183156389 \\ 4 & 3 & -0.0183156389 \\ 4 & 4 & 1 \\ \end{pmatrix}

En rappelant la formule générale des Q_{\alpha},

Q_{\alpha} = \sum_{i = 1} ^n \alpha_i - \frac{1}{2} \cdot \sum_{i_1, i_2} ^n \alpha_{i_1} \alpha_{i_2} Q(i_1, i_2)

, nous obtenons la dérivation suivante,

\alpha_i ^* = \frac{\partial Q_{\alpha}}{\partial \alpha_i} = 1

\Rightarrow \alpha_i ^* = 1 - \sum_{i_2 = 1} ^n \alpha_{i_2} \times Q(i,i_2) = 1

\Rightarrow \alpha_i ^* = - \sum_{i_2 = 1} ^n \alpha_{i_2} \times Q(i,i_2)

Si nous raisonnons en termes de matrice afin de simplifier les calculs, cela revient à transposer les éléments de \mathbf{Q} en fonction de i_1. Nous obtenons alors la matrice,

\begin{pmatrix} \alpha_1 & \alpha_2 & \alpha_3 & \alpha_4 \\ -1 & 0.0183156389 & 0.0183156389 & -0.0003354626 \\ 0.0183156389 & -1 & -0.00032354626 & 0.0183156389 \\ 0.0183156389 & -0.0003354626 & -1 & 0.0183156389 \\ -0.0003354626 & 0.0183156389 & 0.0183156389 & -1 \\ \end{pmatrix}

Comme sous sommes sous contraintes que \sum_{i = 1} ^5 Y_i \times \alpha_i ^* = 0, nous rajoutons la condition au système à résoudre ainsi que la valeur de l’égalité (qui vaut -1 car, pour respecter le passage des équations à la forme matricielle du problème, nous faisons pivoter le coefficient constant associé à la somme des \alpha_i à droite),

\mathbf{S} = \begin{pmatrix} Y_t & \alpha_1 & \alpha_2 & \alpha_3 & \alpha_4 \\ -1 & -1 & 0.0183156389 & 0.0183156389 & -0.0003354626 \\ -1 & 0.0183156389 & -1 & -0.00032354626 & 0.0183156389 \\ -1 & 0.0183156389 & -0.0003354626 & -1 & 0.0183156389 \\ -1 & -0.0003354626 & 0.0183156389 & 0.0183156389 & -1 \\ 0 & 1 & -1 & -1 & 1 \end{pmatrix}

Nous résolvons ce système avec une régression linéaire et obtenons,

\alpha_1 ^* = \alpha_2 ^* = \alpha_3 ^* = \alpha_4 ^* = 1.038

Tous nos coefficients sont différents de 0, nous en déduisons que pour ce noyau les quatre observations jouent un rôle de vecteur support.

Maintenant que les coefficients ont été estimés, nous pouvons déterminer la valeur de prédiction de X_1. Nous déterminons dans un premier temps les coefficients \beta associé aux \alpha ^*,

\beta = (1.038 \times (-1) \times (-1) + 1.038 \times 1 \times (-1) + 1.038 \times 1 \times 1 + 1.038 \times (-1) \times 1, 1.038 \times (-1) \times (-1) + 1.038 \times 1 \times 1 + 1.038 \times 1 \times (-1) + 1.038 \times (-1) \times 1)

= (0,0)

\Rightarrow \beta_0 = Y_1 - X_1 ^t \times \beta = Y_1 = -1

Nous pouvons maintenant déterminer la prédiction de X ^1 sur la base du modèle construit,

f(X_1) = \sum_{i = 1} ^4 \alpha_i ^* Y_i <X_i,X_1> + \beta_0

= 1.038 \times (-1) \times ((-1) ^2 + (-1) ^2) + 1.038 \times 1 \times (-1 \times (-1) + (-1) \times 1) + 1.038 \times 1 \times ((-1) \times 1 + (-1) \times (-1)) + 1.038 \times (-1) \times ((-1) \times 1 + (-1) \times 1) + (-1)

= - 1.038 \times 2 + 1.038 \times 0 + 1.038 \times 0 - 1.038 \times (-2) - 1

= -2.076 + 2.076 - 1

\Rightarrow f(X_1) = - 1

Cas du noyau polynomial:

Soit le noyau polynomial de degré 2, de constante C = 1 et de paramètre \gamma = 1,

K(X_{i_1}, X_{i_2}) = (1 + X_{i_1} ^1 X_{i_2} ^1 + X_{i_1} ^2 X_{i_2} ^2) ^2

Nous calculons pour chaque combinaison d’observations le valeur du noyau,

\mathbf{Q} = \begin{pmatrix} i_1 & i_2 & Y_{i_1} \times Y_{i_2} \times K(X_{i_1}, X_{i_2}) \\ 1 & 1 & 9 \\ 1 & 2 & -1 \\ 1 & 3 & -1 \\ 1 & 4 & 1 \\ 2 & 1 & -1 \\ 2 & 2 & 9 \\ 2 & 3 & 1 \\ 2 & 4 & -1 \\ 3 & 1 & -1 \\ 3 & 2 & 1 \\ 3 & 3 & 9 \\ 3 & 4 & -1 \\ 4 & 1 & 1 \\ 4 & 2 & -1 \\ 4 & 3 & -1 \\ 4 & 4 & 9 \\ \end{pmatrix}

En rappelant la formule générale des Q_{\alpha},

Q_{\alpha} = \sum_{i = 1} ^n \alpha_i - \frac{1}{2} \cdot \sum_{i_1, i_2} ^n \alpha_{i_1} \alpha_{i_2} Q(i_1, i_2)

, nous obtenons la dérivation suivante,

\alpha_i ^* = \frac{\partial Q_{\alpha}}{\partial \alpha_i} = 1

\Rightarrow \alpha_i ^* = 1 - \sum_{i_2 = 1} ^n \alpha_{i_2} \times Q(i,i_2) = 1

\Rightarrow \alpha_i ^* = - \sum_{i_2 = 1} ^n \alpha_{i_2} \times Q(i,i_2)

Si nous raisonnons en termes de matrice afin de simplifier les calculs, cela revient à transposer les éléments de \mathbf{Q} en fonction de i_1. Nous obtenons alors la matrice,

\begin{pmatrix} \alpha_1 & \alpha_2 & \alpha_3 & \alpha_4 \\ -9 & 1 & 1 & -1 \\ 1 & -9 & -1 & 1 \\ 1 & -1 & -9 & 1 \\ -1 & 1 & 1 & -9 \\ \end{pmatrix}

Comme sous sommes sous contraintes que \sum_{i = 1} ^5 Y_i \times \alpha_i ^* = 0, nous rajoutons la condition au système à résoudre ainsi que la valeur de l’égalité (qui vaut -1 car, pour respecter le passage des équations à la forme matricielle du problème, nous faisons pivoter le coefficient constant associé à la somme des \alpha_i à droite),

\mathbf{S} = \begin{pmatrix} Y_t & \alpha_1 & \alpha_2 & \alpha_3 & \alpha_4 \\ -1 & -9 & 1 & 1 & -1 \\ -1 & 1 & -9 & -1 & 1 \\ -1 & 1 & -1 & -9 & 1 \\ -1 & -1 & 1 & 1 & -9 \\ 0 & 1 & -1 & -1 & 1 \end{pmatrix}

Nous résolvons ce système avec une régression linéaire et obtenons,

\alpha_1 ^* = \alpha_2 ^* = \alpha_3 ^* = \alpha_4 ^* = 0.125

Tous nos coefficients sont différents de 0, nous en déduisons que pour ce noyau les quatre observations jouent un rôle de vecteur support.

Maintenant que les coefficients ont été estimés, nous pouvons déterminer la valeur de prédiction de X_1. Nous déterminons dans un premier temps les coefficients \beta associé aux \alpha ^*,

\beta = (0.125 \times (-1) \times (-1) + 0.125 \times 1 \times (-1) + 0.125 \times 1 \times 1 + 0.125 \times (-1) \times 1, 0.125 \times (-1) \times (-1) + 0.125 \times 1 \times 1 + 0.125 \times 1 \times (-1) + 0.125 \times (-1) \times 1)

=(0,0)

\Rightarrow \beta_0 = Y_1 - X_1 ^t \times \beta = Y_1 = -1

Nous pouvons maintenant déterminer la prédiction de X ^1 sur la base du modèle construit,

f(X_1) = \sum_{i = 1} ^4 \alpha_i ^* Y_i <X_i,X_1> + \beta_0

= 0.125 \times (-1) \times ((-1) ^2 + (-1) ^2) + 0.125 \times 1 \times (-1 \times (-1) + (-1) \times 1) + 0.125 \times 1 \times ((-1) \times 1 + (-1) \times (-1)) + 0.125 \times (-1) \times ((-1) \times 1 + (-1) \times 1) + (-1)

= - 0.125 \times 2 + 0.125 \times 0 + 0.125 \times 0 - 0.125 \times (-2) - 1

= -0.250 + 0.250 - 1

\Rightarrow f(X_1) = - 1

\bullet Application informatique:

Macro SAS:

%macro svm(DATA=,Y=,kernel=);

 /* Macro-programme pour le calcul des coefficients associés aux vecteurs supports.
       Les paramètres: – DATA, indique la base de données à traiter,
                       – Y, la variable réponse binaire ou multiclasse,
                       – kernel, le noyau (sigmoidale, radial, polynomial de degré 1-2-3-4).
       Sortie: la table coeffs_&kernel. des coefficients selon le noyau paramétré.

    /* Options nécessaires pour éviter de saturer le log et mettre en évidence les erreurs */
options nonotes spool;

   /* Sortie liste des variables explicatives */
proc contents data = &DATA. (drop =  &Y.) out = biblio noprint;
run;

    /* Listing variables explicatives et nombre */
proc sql noprint;
select distinct(name) into: list_vars separated by  »  » from biblio;
select count(*) into: nb_vars from biblio;
quit;

    /* Sortie effectif total */
proc sql noprint;
select count(*) into: nobs from &DATA.;
quit;

    /* Initialisation de la matrice des Y_i1 Y_i2 K(X_i1,X_i2) */
data Q;
run;

    /* Calcul par paire d’observations */
%do i1 = 1 %to &nobs.;

        /* Récupération 1ère observation */
data obs1;
set &DATA.;
if _n_ = &i1.;
rename &Y. = &Y.1;
%do p = 1 %to &nb_vars.;
rename %scan(&list_vars.,&p.) = %scan(&list_vars.,&p.)_1;
%end;
run;

%do i2 = 1 %to &nobs.;

            /* Récupération 2nde observation */
data obs2;
set &DATA.;
if _n_ = &i2.;
rename &Y. = &Y.2;
%do p = 1 %to &nb_vars.;
rename %scan(&list_vars.,&p.) = %scan(&list_vars.,&p.)_2;
%end;
run;

            /* Calcul de Y_i1 Y_i2 K(X_i1,X_i2) en fonction du noyau paramétré et de la paire d’observations considérées */
data pair;
merge obs1 obs2;
score = 0;
score1 = 0;
score2 = 0;
%do p = 1 %to &nb_vars.;
score = score + %scan(&list_vars.,&p.)_1*%scan(&list_vars.,&p.)_2;
score1 = score1 + %scan(&list_vars.,&p.)_1*%scan(&list_vars.,&p.)_1;
score2 = score2 + %scan(&list_vars.,&p.)_2*%scan(&list_vars.,&p.)_2;
%end;
%if &kernel. = sigmoidale %then %do;
score = &Y.1*&Y.2*tanh(1 + score);
%end;
%if &kernel. = radial %then %do;
score = &Y.1*&Y.2*exp(-sqrt((score1 + score2 – 2 * score)**2));
%end;
%if &kernel. = polynomial1 %then %do;
score = &Y.1*&Y.2*(1 + score);
%end;
%if &kernel. = polynomial2 %then %do;
score = &Y.1*&Y.2*(1 + score)**2;
%end;
%if &kernel. = polynomial3 %then %do;
score = &Y.1*&Y.2*(1 + score)**3;
%end;
%if &kernel. = polynomial4 %then %do;
score = &Y.1*&Y.2*(1 + score)**4;
%end;
run;

data Qadd;
set pair;
Obs1 = &i1.;
Obs2 = &i2.;
keep Obs1 Obs2 score;
run;

data Q;
set Q Qadd;
run;

%end;

%end;

    /* Finalisation de la liste des croisements et de leur score */
data Q;
set Q;
if Obs1 ne .;
score_deriv = – score;
run;

    /* Initialisation de la matrice des lagrangiens par la contrainte à respecter sum Y_i alpha_i = 0 */
proc transpose data = &DATA. (keep = &Y.) out = matsvm;
run;

data matsvm;
set matsvm;
score = 0;
drop _NAME_;
run;

    /* Pour chaque observation, on récupère la dérivée partielle pour le calcul des Lagrangiens */
%do i = 1 %to &nobs.;

data Qt;
set Q;
if Obs1 = &i.;
keep score_deriv;
run;

proc transpose data = Qt out = matsvm_add;
run;

data matsvm_add;
set matsvm_add;
score = -1;
drop _NAME_;
run;

data matsvm;
set matsvm matsvm_add;
run;

%end;

data matsvm;
set matsvm;
run;

    /* Sortie du dernier Lagrangien en guise de borne supérieur des Lagrangiens à considérer */
data _null_;
set &DATA.;
if _n_ = &nobs.;
lim_vars = compress(« COL »||_n_);
call symput (« lim_vars »,lim_vars);
run;

    /* Calcul des Lagrangiens par régression linéaire */
proc reg data = matsvm;
model score = COL1 — &lim_vars. / noint;
ods exclude NObs FitStatistics ParameterEstimates ANOVA;
ods output ParameterEstimates = coeffs;
run;

    /* Finalisation de la table des résultats finaux */
data Coeffs_&kernel.;
set coeffs;
Variable = tranwrd(Variable, »Alpha », »COL »);
Support = « N »;
if Estimate ne 0 then Support = « O »;
keep Variable Estimate Support;
run;

    /* Suppression des tables temporaires inutiles */
proc datasets lib = work nolist;
delete matsvm biblio Q Obs1 Obs2 pair Qadd Qt matsvm_add coeffs;
run;

    /* Reset des options */
options notes nospool;

%mend;

Package et fonction R: https://cran.r-project.org/web/packages/e1071/index.html

\bullet Bibliographie:

– The Elements of Statistical Learning de Trevor Hastie, Robert Tibshirani et Jerome Friedman.

– The Top Ten Algorithms in Data Minig de Xindong Wu et Vipin Kumar.

– Data Mining et statistique décisionnelle. L’intelligence des données de Stéphane Tufféry.

– Probabilité, analyse des données et statistique de Gilbert Saporta.

– SVM, Support Vector Machine. Machines à Vecteurs de Support – Séparateurs à Vaste Marge de Ricco Rakotomalala.

– Méthodes à noyaux – PART I de Alain Rakotomamonjy et Stéphane Canu.

– L’Article web: https://www.math.univ-toulouse.fr/~besse/Wikistat/pdf/st-m-app-svm.pdf

– Sélection de variables par SVM: application à la détection de piétons SVM Variable selection with application to pedestrian detection d’Alain Rakotomamonjy et de Frédéric Suard.

La méthode DISQUAL

add

\bullet Présentation:

La méthode DISQUAL, qui est l’acronyme de DIScrimination sur variables QUALitatives, a été élaboré par Gilbert Saporta en 1975. Elle permet de discriminer une variable polychotomique Y à K \geq 2 classes à partir d’une matrice de variables qualitatives \mathbf{X} = (X ^1, \cdots, X ^P).

Le principe de base de la méthode DISQUAL est de combiner les bienfaits de l’Analyse des Correspondances Mulitples (ACM) avec celles de l’analyse discriminante linéaire ou quadratique de Fisher. Cette caractéristique fait de la méthode DISQUAL une sorte d’équivalent de la régression PLS (qui se base sur un mixte d’Analyse en Composantes Principales et de régression linéaire ou logistique), jouissant ainsi des mêmes propriétés bénéfiques qui sont: la possibilité de conserver plus de variables explicatives dans le modèle, des poids correctement réparties entre ces dernières et l’affranchissement des méfaits de la multicolinéarité.

Si la méthode DISQUAL est à réserver exclusivement à des variables qualitatives, une approche consistant à découper les variables continues en intervalles de valeurs pour les catégoriser est possible afin d’utiliser cet outil sur un format mixte de données.

\bullet DISQUAL:

Hypothèse préliminaire:  Y à K \geq 2 classes et \mathbf{X} qualitative.

L’Algorithme associé à la méthode DISQUAL se déroule à cinq étapes qui sont les suivantes,

– Étape 1: Faire une ACM sur \mathbf{X} et en sortir les projections des variables sur l’ensemble des composantes \mathbf{V} ainsi que celles des observations \mathbf{O}.

– Étape 2: Former la matrice (Y,\mathbf{O}) et sélectionner les axes non nuls et ainsi que ceux séparant au mieux les groupes de Y selon une analyse discriminante de Fisher linéaire ou quadratique sur les axes factorielles. Le critère utilisé en général est le \lambda de Wilks au seuil de 15\% puis 5\% en fonction du nombre d’axes retenus. Suite à ce premier filtre appliqué, (Y,\mathbf{O}) devient (Y,\mathbf{O} ^{D'}), D' ensemble des axes factorielles retenus.

– Étape 3: Appliquer un second filtre sur [latex](Y,\mathbf{O} ^{D'}) en retirant une par une les variables dont la p-valeur est la plus forte et en usant du critère de l'Aire (AUC), du volume (VUS) ou hypervolume (HUM) sous la Courbe ROC (AUC) en fonction du nombre de classe de Y obtenu selon une analyse discriminante de Fisher linéaire ou quadratique sur les axes factorielles. Nous supprimons ainsi les dimensions à partir desquelles ce critère chute. Suite à ce second filtre appliqué, (Y,\mathbf{O} ^{D'}) devient (Y,\mathbf{O} ^D)

- Étape 4: Nous lançons une analyse discriminante de Fisher linéaire ou quadratique sur (Y,\mathbf{O} ^D) et obtenons ainsi les coefficients associés à chacune de nos composantes. Nous recherchons l'estimation combinée à partir de ces coefficients obtenus en soustrayant les coefficients par composantes en fonction des groupes de Y (dans le cas où Y est polychotomique, la distance euclidienne peut être une option).

- Étape 5: Le calcul des coefficients \beta du modèle prédictif se fait alors en récupérant \mathbf{V} ^D et en multipliant les coordonnées de nos modalités de variables par les coefficients respectifs à chacun de ces axes et obtenus en étape 4. La règle décisionnelle est alors de la forme suivante,

Y = \sum_{p = 1} ^P \sum_{m_p = 1} ^{\sharp \lbrace \mbox{modalites de} X ^p \rbrace} \beta_p 1|_{X ^p = m_p}

A noter que nous pouvons également fixer une étape de sélection des variables à cet algorithme lors de l'étape 1. Il faudrait alors lancer une première ACM sur \mathbf{X}, supprimer les variables dont l'ensemble des modalités ne sont que très peu contributives à la formation des axes et relancer une nouvelle ACM sur l'ensemble épuré de ces variables.

\bullet Annexe théorique:

Nous présentons dans cette partie de l'article une esquisse de la démarche méthodologique justifiant l'utilisation de la méthode DISQUAL.

La méthode DISQUAL se base sur un mélange d'ACM et d'analyse discriminante de Fisher. L'Utilisation d'une ACM sur un tableau de données qualitative \mathbf{X} revient à, dans un premier temps, transformer notre matrice en un tableau disjonctif complet, soit un tableau de Q indicatrices valant 1 si l'observation a choisi tel modalité de tel variable et 0 sinon.

Le problème de discrimination revient alors à lancer un analyse discriminante de Fisher sur ce tableau disjonctif et à estimer les coefficients propres à chaque modalité de chaque variable (ce qui revient à dire à chaque variable du dit-tableau). Nous pouvons alors reformuler l'objectif en la maximisation de la distance de Mahalanobis entre les centres de gravité des différents groupes de Y.

En ce sens, la méthode DISQUAL s’inspire donc très largement de la régression sur composantes principales (ou de la régression PLS) en lançant une analyse discriminante de Fisher sur nos axes factoriels.

Soit b ^q les composantes de l'ACM et \lambda_q les valeurs propres qui leur sont associés. L'Une des formidables propriétés des axes factoriels construis via une ACM est qu'ils sont orthogonaux, rendant l'inversion des matrices nécessaires au calcul possible comparé au cas où nous aurions directement travaillé sur le tableau disjonctif complet. Nous pouvons alors obtenir le score \delta issue de la fonction de Fisher (dans le cas où Y binaire),

\delta = \sum_{q = 1} ^Q u_q b ^q

, avec,

\mathbf{u} = \begin{pmatrix} \cdot \\ u_q \\ \cdot \\ \end{pmatrix} = \mathbf{V} ^{-1} (\mu_1 - \mu_2) = \begin{pmatrix} \cdot \\ \frac{\overline{b_1} ^q - \overline{b_2} ^q}{\lambda_q} \\ \cdot \\ \end{pmatrix}

L'Avantage de s'appuyer sur l'ACM est qu'il suffit d’effectuer la combinaison linéaire avec les mêmes coefficients u_q des coordonnées de ses catégories. En effet, nous avons b ^q = \mathbf{X} a ^q, où a ^q vecteur des coordonnées des modalités dans le plan factoriel. Nous avons alors,

\delta = \sum_{q = 1} ^Q u_q \mathbf{X} a ^q = \mathbf{X} \sum_{q = 1} ^Q u_q a ^q

Le score s'exprime alors comme la somme de la combinaison linéaire des indicatrices des modalités pour chaque individu. A noter que la règle décisionnelle ne comporte volontairement pas de terme constant.

\bullet Exemple:

Soit le jeu de données ci-dessous,

add

Appliquons la méthode DISQUAL, nous commençons donc par une ACM sur \mathbf{X} = (X ^1, X ^2, X ^3, X ^4, X ^5) et sortons les coordonnées de nos individus sur les axes factoriels conçus,

add

Nous pouvons déjà filtrer six composantes sur les quinze attendues puisque nulles. Nous conservons donc Dim1, Dim2, \cdots, Dim9.

Maintenant nous devons filtrer les dimensions à partir du modèle Y = \beta_0 + \sum_{d = 1} ^D \beta_d Dim_d. Nous utilisons pour cela le test du \lambda de Wilks et obtenons les résultats suivants classés par ordre d'intégration selon une méthode descendante,

\begin{tabular} {|l|c|} \hline Axe & p-valeur \\ \hline Dim1 & 0.0007 \\ \hline Dim2 & 0.0010 \\ \hline Dim4 & 0.0078 \\ \hline Dim3 & 0.0415 \\ \hline Dim9 & 0.0372 \\ \hline Dim5 & 0.0198 \\ \hline Dim8 & 0.0284 \\ \hline \end{tabular}

Au seuil de 5 \%, sept axes sont significatifs. Par conséquent, nous les conservons tous à cette étape. Appliquons maintenant un second filtre basé sur l'AUC. Nous obtenons alors,

\begin{tabular} {|l|c|} \hline Suppression & AUC \\ \hline Dim8 & 1 \\ \hline Dim5 & 1 \\ \hline Dim9 & 1 \\ \hline Dim3 & 1 \\ \hline Dim4 & 0.99 \\ \hline Dim2 & 0.92 \\ \hline Dim1 & 0.51 \\ \hline \end{tabular}

Il ressort du tableau ci-dessus que nous pouvons supprimer les dimensions Dim_8, Dim_5, Dim_9, Dim_3) puisque sur le modèle Dim_4, Dim_2, Dim_1, l'AUC reste maximale.

Maintenant nous faisons tourner une nouvelle analyse discriminante de Fisher sur le modèle retenu et calculons le score en fonction des prédictions. Nous obtenons,

\begin{tabular} {|l|c|c|c|} \hline Axe & Coeff A & Coeff B & Score \\ \hline Cst & -2.18810 & -2.18810 & 0 \\ \hline Dim1 & 4.77941 & -4.77941 & 9.5588283314 \\ \hline Dim2 & 3.87168 & -3.87168 & 7.7433688952 \\ \hline Dim4 & -3.53072 & 3.53072 & -7.061434395 \\ \hline \end{tabular}

Nous avons maintenant nos coefficients par axe factoriel retenu. Place à l'estimation des coefficients de la règle décisionnelle en fonction des coordonnées des modalités de variables sur les axes Dim1, Dim2, Dim4,

\begin{tabular} {|l|c|c|c|c|} \hline Modalites & Dim1 & Dim2 & Dim4 & Coefficient \\ \hline X1 = 1m & 1.4525 & -0.0724 & -0.0991 & 14.024033312 \\ \hline X1 = 2m & -0.5044 & 1.0807 & -0.4710 & 2.507363696 \\ \hline X1 = 3m & -0.7407 & -1.0187 & -0.0623 & -14.52796362 \\ \hline X2 = 1m & 0.0672 & 0.1352 & 1.1090 & -6.141877623 \\ \hline X2 = 2m & 0.4744 & -0.59144 & -0.2168 & 1.4867091678 \\ \hline X2 = 3m & -0.5321 & 0.4756 & -0.7337 & 3.7777573659 \\ \hline X3 = 1m & -0.7178 & 0.4918 & 0.7324 & -8.225486258 \\ \hline X3 = 2m & 1.0569 & 0.4305 & -0.0022 & 13.451550502 \\ \hline X3 = 3m & -0.4416 & -0.8520 & -0.6256 & -6.401133709 \\ \hline X4 = 1m & 1.4525 & -0.0724 & -0.0991 & 14.024033312 \\ \hline X4 = 2m & -0.6846 & 0.6239 & -0.4711 & 1.163251329 \\ \hline X4 = 3m & -0.5604 & 0.5619 & 0.5560 & -13.63385131 \\ \hline X5 = 1m & -0.2261 & 1.2179 & -0.2590 & 9.098504086 \\ \hline X2 = 2m & 0.7383 & -0.1965 & 0.1248 & 4.6541658144 \\ \hline X5 = 3m & -0.5445 & 0.8474 & 0.0972 & -12.4528836 \\ \hline \end{tabular}

\bullet Application informatique:

Procédure SAS: http://od-datamining.com/download/#macro

Package et fonction R: http://finzi.psych.upenn.edu/library/DiscriMiner/html/disqual.html

\bullet Bibliographie:

- Probabilité, analyse des données et Statistique de Gilbert Saporta.

- Data Mining et statistique décisionnelle de Stéphane Tufféry.

Les méthodes de boosting

add.png

\bullet Présentation:

Le boosting est une famille d’outils d’analyse discriminante très puissants. L’idée est d’améliorer la règle de discrimination d’un classifieur sur plusieurs sous-parties de l’hyperplan formé par (Y, X ^1, \cdots, X ^P) et d’établir la classification finale au travers d’un vote global pondéré. Chaque partie du plan se voit alors attribuer une pondération en fonction de sa performance de discrimination et qui est mise à jour itération par itération afin d’adapter la frontière décisionnelle en fonction des zones les plus difficiles à classer.

Le classifieur utilisé est la plupart du temps la régression linéaire, logistique ou les arbres de décision. Ils portent le nom de « classifieurs fiables » à la fois car ils sont simples d’utilisation mais également car l’objectif du boosting est d’améliorer la règle de discrimination lorsqu’elle est peu performante.

Le boosting appartient à la catégorie des méthodes d’apprentissage statistique nécessitant un jeu de d’apprentissage pour la conception de la règle d’apprentissage et un jeu de test pour sa validation afin de limiter les risques de sur-apprentissage. Comme la plupart d’entre elles, le boosting a le principal défaut de concevoir une règle décisionnelle peu lisible d’où son surnom de boîte noire.

Parmi les méthodes de boosting, la plus célèbre reste l’algorithme AdaBoost conçu par Yoav Freun et Robert Shapire en 1997 et définit pour la discrimination d’une variable réponse Y binaire à partir d’un ensemble de variables continues explicatives \mathbf{X} = (X ^1, \cdots, X ^P).

D’autres algorithmes de boosting sont particulièrement répandus, nous pouvons ainsi citer: l’algorithme Gentle Adaboost, Real Adaboost, MART et RankBoost.

\bullet Les algorithmes de Boosting:

Hypothèse préliminaire: Y binaire ou continue, \mathbf{X} continue.

L’algorithme générale:

Avant d’aborder les algorithmes Adaboost et MART, nous présentons l’algorithme de boosting de base.

– Étape initiale: définir un classifieur faible f et un nombre T d’itérations d’entrainement pour l’optimisation des pondérations.

– Étape 1: Soit \mathbf{E}_1 = \mathbf{E} = [Y \mathbf{X}], la matrice de données d’origine et pour laquelle chaque observation a un poids similaire aux autres et égal à 1.

– Étape 2: Pour t \in \lbrace 2, \cdots, T \rbrace,

—- Nous appliquons notre classifieur faible: h_t = f(\mathbf{E}_1)

—- Nous mesurons le taux d’erreurs de classification de la règle décisionnelle conçue précédemment: \epsilon_t = Pr_{x ~ \mathbf{E}_t, y} I[h_t(x) \neq y]

—- Nous mettons à jour les pondérations et chaque observation mal classée se voit attribuer un poids plus important. Ainsi, \mathbf{E}_t devient \mathbf{E}_{t + 1}

– Étape 3: Nous déterminons la règle décisionnelle finale: H(x) = \mbox{Vote a partir des } h_t(x), \forall t \in [1,T]

L’algorithme AdaBoost

L’Algorithme AdaBoost a été pensé pour le cas où Y est binaire. Il se décline de la manière suivante:

– Étape initiale: initialisation des pondérations associées aux observations,

w_i ^1 = \frac{1}{N}, \forall i \in [1,N]

– Étape 1: \forall t \in [2,T]

—- Application du classifieur paramétré sur le jeu d’apprentissage en utilisant les pondérations w_i ^{t - 1}, \forall i \in [1,N], nous obtenons la règle décisionnelle f_t(x).

—- Calcul du taux d’erreurs de la règle décisionnelle définie précédemment,

\epsilon_t = Pr_{x ~ \mathbf{E}_t, y} I[f_t(x) \neq y]

—-Si \epsilon_t >0.5 alors le classifieur fait mieux qu’une règle aléatoire et nous conservons le modèle, sinon nous ajustons,

  • le coefficient de pondération:

\alpha_t = \frac{1}{2} ln(\frac{1 - \epsilon_t}{\epsilon_t})

  • les pondérations:

w_i ^t = w_i ^{t - 1} e ^{\alpha_t \cdot I_{(y_i \neq f_t (x_i))}}, \forall i \in [1,N]

– Étape 2: Nous obtenons le modèle final,

H(x) = sign[\sum_{t = 1} ^T \alpha_t f_t (x)]

L’algorithme Multiple Additive Regression Trees (MART)

L’Algorithme MART se base sur l’utilisation d’arbres de régression. Par conséquent, il est adapté au cas où Y est continue. Les différentes étapes d’estimation sont les suivantes:

– Étape initiale: Initialiser f_0 (x) = arg min_{\delta} \sum_{i = 1} ^N L(y_i, \delta), avec L fonction de perte à minimiser. Et détermination du nombre d’itérations T.

– Étape 1: Pour t \in [1, T],

— Pour i \in [1, N] calculer,

r_{i,t} = - [\frac{\partial L(y_i, f(x_i))}{\partial f(x_i)}]_{f = f_{t - 1}}

— Estimer le modèle par arbre de régression sur r_{i,t} et donnant les régions décisionnelles R_{j,t}, \forall j \in \lbrace 1, \cdots, J_t \rbrace

— Pour j \in \lbrace 1, \cdots, J_t \rbrace, calculer,

\delta_{j,t} = arg min_{\delta} \sum_{x_i \in R_{j,t}} L(y_i, f_{t - 1} (x_i) + \delta)

— Mettre à jour f_t (x) = f_{t - 1} (x) + \sum_{j = 1} ^{J_t} \delta_{j,t} I (x \in R_{j,t})

– Étape 2: Nous avons enfin la règle décisionnelle finale de forme: \hat{f} (x) = f_T (x)

\bullet Annexe théorique:

Cette partie de l’article présente une esquisse de la justification de l’usage de la fonction de perte exponentielle.

L’algorithme AdaBoost se base sur l’idée nouvelle d’une fonction de perte exponentielle dont le principal avantage est computationnel. Ainsi, nous montrons facilement que,

f ^* (x) = arg min_{f(x)} E_{Y|x} (e ^{-Y f(x)}) = \frac{1}{2} log \frac{P(Y = 1 | x)}{P(Y = -1 | x)}

Et équivalent à,

P(Y = 1 | x) = \frac{1}{1 + e ^{-2 f^* (x)}}

De cette manière, l’expansion additive produite par l’algorithme AdaBoost revient à estimer une partie du log-odds de P(Y = 1 | x). Ce qui justifie l’usage de son signe comme règle de classification décrite par la formule,

G(x) = sign (\sum_{m = 1} ^M \alpha_m G_m (x))

Un autre critère de minimisation de perte est la log-vraisemblance de la loi binomial négative (appelée également déviance ou entropie croisée), interprétant f comme une transformation logit. Nous avons,

p(x) = P(Y = 1 | x) = \frac{e ^{f(x)}}{e ^{-f(x)} + e ^{f(x)}}  = \frac{1}{1 + e ^{-2 f(x)}}

, avec Y ' = \frac{Y + 1}{2} \in \lbrace 0, 1 \rbrace. La fonction de perte est alors,

l(Y,p(x)) = Y ' log (p(x)) + (1 - Y ') log(1 - p(x))

Et de déviance,

-l (Y, f(x)) = log (1 + e ^{-2Y f(x)})

De la formule précédente, nous voyons que,

E_{Y|x} [-l(Y,f(x))] = E_{Y|x} [e ^{-Y f(x)}]

\bullet Application informatique:

Procédure SAS: ND.

Package et fonction Rhttps://www.rdocumentation.org/packages/JOUSBoost/versions/2.1.0/topics/adaboost

\bullet Bibliographie:

– The Elements of Statistical Learning de Trevor Hastie, Robert Tibshirani et Jerome Friedman.

– The top ten algorithms in data mining de Xindong Wu et Vipin Kumar.

– Probabilité, analyse des données et Statistique de Gilbert Saporta.

– Data mining et statistique décisionnelle. L’intelligence des données de Stéphane Tufféry.

Le modèle additif généralisé

add.png

\bullet Présentation:

Le modèle additif généralisé, que nous pouvons retrouvons sous l’acronyme GAM, a été élaboré par Trevor Hastie et Rob Tibshirani en 1990.  Le modèle GAM se base sur la somme (doù le terme d’additif) de fonctions f_1, \cdots, f_P de transformation des différentes variables explicatives \mathbf{X} = (X ^1, \cdots, X^P).

La variable à expliquer Y peut être de tout format et les variables explicatives \mathbf{X} = (X ^1, \cdots, X ^P) doivent êtres continues.

Le modèle GAM nécessite un jeu d’apprentissage pour sa construction et un jeu de test pour sa validation car il est sujet à un possible fort sur-apprentissage.  En dépit de sa puissance avérée pour modéliser un phénomène, il reste peu conseillé si nous disposons d’un nombre P de variables explicatives trop important du fait d’un coût de calcul très élevé.

\bullet Le modèle additif généralisé:

Hypothèse préliminaire: \mathbf{X} continue.

Le modèle:

La forme générale du modèle est,

E[Y | X ^1, \dots, X ^P] = \alpha + \sum_{p = 1} ^P f_p (X ^p)

Avec f_p, \forall p \in [1, P], fonction de transformation associée à la variable X ^p ou au uplet de variables explicatives. Deux principaux algorithmes d’estimation des paramètres sont utilisés: l’algorithme Backfitting et l’algorithme général local de scoring,

L’Algorithme backfitting (Friedman et Stuezle, 1981):

– Etape d’initialisation: poser f_0 = E[Y] = \frac{1}{n} \sum_{i = 1} ^n Y_i, et f_p ^0 = 0, \forall p \in[1, P]

– Etape itérative t: pour p \in [1,P], déterminer,

  • R_p = Y - f_0 - \sum_{k = 1} ^{p - 1} f_k ^t (X ^k) - \sum_{k = p + 1} ^P f_k ^{t - 1} (X ^k)
  • f_p ^t = E[R_p | X ^p]

Tant que,

\frac{1}{n} || Y - s_0 - \sum_{p = 1} ^P s_p ^t (X ^p) || ^2 décroît ou satisfait le critère de convergence

A noter que pour p = 1, \sum_{k = 1} ^{1 - 1} s_k ^t (X ^k) = \sum_{k = 1} ^0 s_k ^t (X ^k) = 0 et pour p = P, \sum_{k = P + 1} ^P s_k ^t (X ^k) = 0.

De plus, l’estimation des différentes fonctions de lissage à l’instant t, f_p ^t, se fait selon la fonction lien retenue (modèle linéaire, logistique, etc).

Les transformations:

Le modèle GAM se base donc sur la transformation de variables ou de combinaisons de variables afin d’optimiser le modèle prédictif construit. A ce titre, trois principales fonctions existent: La fonction Cubic Smoothing Spline, la fonction de régression locale (LOESS) et la fonction Thin-Plate Smoothing Spline.

– La fonction Cubic Smoothing Spline de paramètre \lambda \in [0, + \infty[:

RSS(f,\lambda) = \sum_{i = 1} ^n [y_i - f(x_i)] ^2 + \lambda \int (f '' (t)) ^2 dt

– La fonction LOESS de paramètre le noyau K:

min_{\alpha(x_0), \beta(x_0)} \sum_{i = 1} ^n K_{\lambda} (x_0,x_i) [y_i - \alpha(x_0) - \beta(x_0) x_i] ^2

– La fonction Thin-Plate Smooting Spline:

f(x) = \beta_0 + \beta ^t x + \sum_{i = 1} ^n \alpha_i h_i (x), h_i (x) = \theta(|| x - x_p||), \theta(z) = z ^2 log z ^2

\bullet Annexe théorique:

Cette partie de l’article présente une esquisse de la démonstration des fonctions à estimer pour déterminer les valeurs de \mathbf{X} selon la projection f = (f_1, f_2, \cdots).

– Nous avons vu que la fonction Cubic Smoothing Spline de paramètre \lambda \in [0, + \infty[ était définie ainsi:

RSS(f,\lambda) = \sum_{i = 1} ^n [y_i - f(x_i)] ^2 + \lambda \int (f '' (t)) ^2 dt

Lorsque \lambda = 0 alors f peut être n’importe quel fonction d’interpolation de X. Si \lambda = \infty, alors f est la fonction des moindres carrés partiels et aucune derivée seconde ne peut être obtenue.

Pour obtenir une version simplifiée de la fonction Cubic Smoothing Spline, nous posons D_j(x) un sous-ensemble de D fonctions basiques et \theta un vecteur de paramètres à estimer. La fonction d’origine peut alors se réécrire:

RSS(\theta,\lambda) = (y - D \theta) ^t (y - D \theta) + \lambda \theta ^t \Omega_N \theta

, avec \lbrace D \rbrace_{i,j} = D_j (x_i) et \lbrace \Omega_N \rbrace_{j,k} = \int D_j '' (t) D_k '' (t) dt. La solution est alors,

\hat{\theta} = (D ^t D + \lambda \Omega_D) ^{-1} D ^t y

, qui s’apparente à la régression ridge généralisée. Après développement, nous obtenons ainsi la forme simplifiée de f,

\hat{f}(x) = \sum_{j = 1} ^D D_j (x) \hat{\theta}_j

– Nous avons vu que la fonction LOESS de paramètre le noyau K était définie ainsi:

min_{\alpha(x_0), \beta(x_0)} \sum_{i = 1} ^n K_{\lambda} (x_0,x_i) [y_i - \alpha(x_0) - \beta(x_0) x_i] ^2

La fonction d’estimation est de la forme,

\hat{f} (x_0) = \hat{\alpha} (x_0) + \hat{\beta} (x_0) x_0

Nous définissons alors la fonction du vecteur de valeurs: b(x) ^t = (1,x), \mathbf{B} la matrice de régression de taille N \times 2 composée des b(x_i) ^t et \mathbf{W}(x_0) la matrice diagonale de taille N \times N et dont l’élément i est K_{\lambda} (x_0, x_i). Nous avons alors,

\hat{f} (x_0) = b(x_0) ^t \mathbf{B} ^t \mathbf{W} (x_0) \mathbf{B} ^{-1} \mathbf{B} ^t \mathbf{W} (x_0) y = \sum_{i = 1} ^N l_i (x_0) y_i

Ce qui nous donne une expression explicite pour l’estimation de la LOESS.

– Nous avons vu que la fonction Thin-Plate Smooting Spline était définie ainsi:

f(x) = \beta_0 + \beta ^t x + \sum_{i = 1} ^n \alpha_i h_i (x)

, avec h_i (x) = \theta(|| x - x_p||) et \theta(z) = z ^2 log z ^2.

Soit \mathbf{X} \in R ^2 et les fonctions de base de la forme h_{1,k} (X ^1), \forall k \in 1, \cdots, M_1 des coordonnées de X ^1 et de même pour la fonction h_{2,k} (X ^2), \forall 1, \cdots, M_2. Nous définissons alors,

g_{j,k} (\mathbf{X}) = h_{1,j} (X ^1) h_{2,k} (X ^2), j = 1, \cdots, M_1 et k = 1, \cdots, M_2

, que nous pouvons généraliser avec la forme suivante,

g(\mathbf{X}) = \sum_{j = 1} ^{M_1} \sum_{k = 1} ^{M_2} \theta_{j,k} g_{j,k} (\mathbf{X})

Pour estimer cette fonction, l’idée est de résoudre le problème suivant,

min_J \sum_{i = 1} ^N \lbrace y_i - f(x_i) \rbrace ^2 + \lambda J[f]

, avec J fonction de pénalisation visant à stabiliser f.

De plus, si le paramètre \lambda \rightarrow 0 alors la solution s’approche d’une fonction d’interpolation et si \lambda \rightarrow \infty alors elle s’approche d’un plan des moindres carrés. Pour les valeurs intermédiaires de \lambda, la solution peut être représentée comme une fonction linéaire dont les coefficients sont obtenus par régression ridge.

La solution est alors de la forme,

f(x) = \beta_0 + \beta^t x + \sum_{j = 1} ^N \alpha_j h_j (x)

, où h_j(x) = \eta(||x - x_j||) et \eta(z) = z ^2 log z ^2.

\bullet Application informatique:

Procédure SAS: https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_gam_sect004.htm

Package et fonction R: https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/gam.html

\bullet Bibliographie:

– The Elements of Statistical Learning de Trevor Hastie, Robert Tibshirani et Jerome Friedman.

– Data Mining et statistique décisionnelle. L’Intelligence des données de Stéphane Tufféry.

https://www.math.univ-toulouse.fr/~besse/pub/Appren_stat.pdf

La régression linéaire

add.png

\bullet Présentation:

La régression linéaire permet de modéliser une variable réponse Y continue à partir d’une matrice de P variables explicatives \mathbf{X} = (X ^1, \cdots, X ^P) continues. Une extension de la régression linéaire permet de traiter le cas où \mathbf{X} est exclusivement formée de variables qualitatives, nous parlons alors de régression linéaire sur variables modales.

Difficile de retrouver le véritable inventeur de la régression linéaire, de plus il est compliqué de distinguer celui qui en a trouver la forme théorique générale de celui qui a inventé l’une des trois méthodes d’estimation des paramètres utilisées. Pour faire un rapide historique, la littérature semble s’accorder sur le fait que c’est à Roger Joseph Boscovich qui en est le précurseur en 1757. Ses travaux seront repris par Pierre Simon de Laplace en 1789, Carl Friedrich Gauss en 1795 et Adrien Marie Legendre qui en publiera la méthode des moindres carrés en 1805.

La régression linéaire se base sur une modélisation de type linéaire et demeure la méthode de prédiction la plus utilisée étant donné sa simplicité de mise en œuvre. Cinq conditions sont à tester au préalable pour garantir une utilisation optimale de la régression linéaire:

– l’absence de multi-colinéarité au sein de \mathbf{X},

– l’indépendance des résidus \epsilon_i, \forall i \in [1,n], entre eux,

– l’exogénéité des variables explicatives X ^p, \forall p \in [1,P] soit l’indépendance entre elles et les résidus \epsilon_i, \forall i \in [1,n],

– l’homoscédasticité des résidus \epsilon soit que leur variance est constante,

– la normalité des résidus \epsilon,

A l’exception de la première hypothèse, et dans son cas extrême, l’utilisation de la régression linéaire est possible, sauf que plus nous nous éloignons des cinq hypothèses et plus nous risquons d’obtenir des estimateurs biaisés et donc peu efficaces.

Enfin, trois principales extensions à la régression linéaire existent: la régression RIDGE, la régression LASSO et la régression ELASTIC NET visant à améliorer l’estimation des paramètres en fonction des spécificités de \mathbf{X}. Ces outils portent également le nom commun de régression pénalisée (étant donné que l’idée est de se soustraire à la grande variance des estimateurs en positionnant un paramètre de pénalisation visant à réduire cette variabilité).

\bullet La régression linéaire:

Hypothèses préliminaires: Variable réponse Y continue. Variables explicatives \mathbf{X} continues. Pas de multicolinéarité et exogénéité de \mathbf{X}; et indépendance, homoscédasticité et normalité des résidus.

Le modèle:

Le modèle associé à la régression linéaire s’exprime au travers de la formule suivante,

Y = \beta_0 + \beta_1 X ^1 + \cdots + \beta ^P X ^P + \epsilon

L’Idée est donc de reconstruire Y à partir des X ^p, p \in [1, P] par une formule ou lien linéaire.

L’estimation des coefficients:

Il existe trois principales méthodes d’estimations des coefficients: les moindres carrés partiels, le maximum de vraisemblance et l’inférence bayésienne. La plus populaire de toute demeure celle des moindres carrés qui permet de fournir des estimateurs \beta non biaisés.

L’estimation des coefficients de régression \beta = (\beta_0, \beta_1, \cdots, \beta_P) se fait alors au travers du calcul,

\beta = (\mathbf{X}_1 ^t \mathbf{X}_1) ^{-1} \mathbf{X}_1 ^t Y

, où \mathbf{X}_1 = (1, \mathbf{X}) la matrice des variables explicatives à laquelle nous rajouter en première colonne le vecteur unitaire dans le but de pouvoir estimer le coefficient constant \beta_0.

Les indicateurs de performance de la régression linéaire:

En notant,

\epsilon = (\epsilon_1, \cdots, \epsilon_n) = (Y_1 - \beta_0 - \sum_{p = 1} ^P \beta_p X_1 ^p, \cdots, Y_n - \beta_0 - \sum_{p = 1} ^P \beta_p X_n ^p)

, autrement dit la différence entre valeurs réelles de Y et valeurs prédites \hat{Y} par le modèle construit via régression linéaire, nous pouvons désigner trois principaux indicateurs régulièrement retrouvés dans la littérature: la vraisemblance, le R ^2 et le R ^2 ajusté noté [atex]R_{adj} ^2[/latex].

  • La vraisemblance,

L = L(\beta_0, \cdots, \beta_P, Y_1, \cdots, Y_n) = \frac{1}{\sqrt{2 \pi \sigma ^2}} e^{-\frac{1}{2 \sigma ^2} \sum_{i = 1} ^n \epsilon_i ^2}

, plus la vraisemblance est forte, meilleur est le modèle.

Le R ^2 (ou coefficient de détermination) offre un indicateur statistique de la modélisation du modèle et notamment de sa capacité à coller aux données réelles utilisées pour le construire. Le R ^2 varie entre [0,1] et plus il s’approche de sa borne supérieure, meilleur est le modèle. Il existe toute une flopée d’indicateurs de type R ^2, en voici les principaux.

  • Le R ^2, de formule,

R ^2 = 1 - \frac{\sum_{i = 1} ^n \epsilon_i ^2}{SST}

  • Le R_{adj} ^2de formule,

R_{adj} ^2 = 1 - \frac{(n - \sharp \lbrace \beta_0 \rbrace) (1 - R ^2)}{n - P - 1}

Avec, pour les deux derniers indicateurs présentés,

SST = \sum_{i = 1} ^n \epsilon_i ^2 + \sum_{i = 1} ^n (\hat{Y}_i - \overline{Y}) ^2

, si le modèle contient une constante (\sharp \lbrace \beta_0 \rbrace = 1). Si le modèle n’en contient pas de constante (\sharp \lbrace \beta_0 \rbrace = 0), alors,

SST = \sum_{i = 1} ^n \epsilon_i ^2 + \sum_{i = 1} ^n \hat{Y}_i ^2

D’autres versions du R ^2 existent comme celui de Nagelkerke, d’Efron, de Cox-Snell, McKelvey-Zavoina, de McFadden et sa version ajustée ou les versions empiriques « count » et « adjusted count ». Cependant ils restent rarement utilisés.

La sélection de variables:

Nous dénombrons trois principaux indicateurs permettant de sélectionner les variables les plus intéressantes dans l’objectif d’optimiser un modèle construit par régression linéaire. Il s’agit du Critère d’Information d’Akaike, noté AIC, le Critère d’Information Bayésien, noté BIC, et le critère Combinatoire des P variables de Mallows, noté C_p de Mallows.

L’utilisation de ces trois critères se fait lors des célèbres algorithmes sélectifs ascendante (FORWARD), descendante (BACKWARD) et mixte (STEPWISE) qui consiste à,

– FORWARD: Partir, pour l’étape initiale, du modèle constant et, en mesurant l’apport des différentes variables rajoutées l’une après l’autre, retenir celle qui apporte le plus d’information et l’intégrer au modèle. Les étapes suivantes procèdent de même façon sauf qu’elle compare le modèle à l’itération it et l’apport potentiel des variables non intégrées à celui-ci. L’Algorithme s’arrête si plus aucune variable n’apporte significativement pas à la qualité du modèle en cours.

– BACKWARD: Partir, pour l’étape initiale, du modèle complet (soit avec les P variables) et, en retirant une à une les P-1 variables, retenir celle qui enlèvle le moins d’information et la retirer du modèle. Les étapes suivantes procèdent de même façon sauf qu’elle compare le modèle à l’itération it et le nerf des variables restantes à celui-ci. L’Algorithme s’arrête si toutes les variables restantes diminuent significativement la qualité du modèle en cours.

– STEPWISE: il s’agit d’un mixte entre l’algorithme FOWARD et l’algorithme BACKWARD qui part du modèle constant, applique une étape FORWARD et enchaîne sur une étape BACKWARD avant de valider le modèle et poursuivre ses itérations. Ainsi, si la variable testée apporte de manière significative au modèle et si le fait de la retirer diminue trop son pouvoir prédictif, alors elle est conservée.

Présentons maintenant les trois critères avec \sigma ^2 = var(\epsilon) est la variance de l’erreur de prédiction,

– Le critère AIC, inventé par Hirotugu Akaike en 1974 et de formule,

AIC = n ln(\frac{\sum_{i = 1} ^n \epsilon_i ^2}{n}) + 2 (P + 1)

, plus la valeur d’AIC est faible et plus le modèle est de bonne qualité.

– Le critère BIC, inventé par Gideon Schwarz en 1978 et de formule,

BIC = n ln(\frac{\sum_{i = 1} ^n \epsilon_i ^2}{n}) + 2 (P + 3) \frac{n \sigma ^2}{SSE} - 2 (\frac{n \sigma ^2}{SSE}) ^2

, plus la valeur du BIC est faible et plus le modèle est de bonne qualité.

– Le C_p de Mallows, inventé Collin Lingwood Mallow et de formule,

C_p = \frac{\sum_{i = 1} ^n \epsilon_i ^2}{\sigma ^2} - n + 2 P

, plus la valeur du C_p de Mallows est faible et plus le modèle est de bonne qualité.

La régression RIDGE:

La régression RIDGE, ou L_2-pénalisation, a été inventée par A. E. Horel en 1962 et consiste à améliorer l’estimation des coefficients \beta lorsque nous sommes en présence de forte multi-colinéarité pouvant apporter un biais substantielle.

L’idée de la régression RIDGE est de minimiser l’erreur quadratique moyenne des estimateurs  au travers de la recherche d’un compromis biais/variance. La méthode permet alors de calculer et d’utiliser des estimateurs biaisés mais de variance plus faible que ceux estimés par méthode des moindres carrés.

L’estimation des coefficients de régression associés au modèle RIDGE se fait alors par la formule,

\beta_R = (\mathbf{X_1} ^t \mathbf{X_1} + \lambda I) ^{-1} \mathbf{X_1} ^t Y

, avec I matrice identité de taille P \times P et lambda paramètre de L_2-pénalisation.

add.png

L’estimation du paramètre de pénalisation \lambda se fait au travers de la formule,

\lambda = \frac{P \sigma ^2}{\beta^t \beta}

La procédure, inspirée par A. E. Hoerl et R. W Kennard en 1970, pour déterminer les paramètres \sigma, \gamma est,

– Lancer l’estimation du modèle en posant \lambda = 0 afin d’obtenir une première estimation du vecteur \beta_R = \beta puisque cela revient à procéder à une régression linéaire standart.

– Calculer l’écart-type des résidus \epsilon_i, i \in [1,n], soit la différence entre Y et les prédictions \hat{Y} via le modèle construit pour \lambda = 0.

Une fois ces deux paramètres estimés, nous pouvons calculer \lambda et en relançant le modèle nous obtenons les coefficients de régression RIDGE finaux \beta_R.

Une autre méthode d’estimation du paramètre \lambda existe, elle consiste à définir un intervalle de variation et à traçer les valeurs des différents coefficients. La valeur minimale pour laquelle les coefficients convergent est la valeur \lambda à retenir.

La régression RIDGE présente l’avantage d’accroître le pouvoir explicatif du modèle en présence de fortes multicolinéarité en attribuant aux variables très informatives et fortement corrélées le même coefficient. Mais aussi le désavantage de ne pas permettre de distinguer les variables les plus contributives des variables les moins contributives à la qualité prédictive du modèle.

A noter qu’en fonction de la littérature, \mathbf{X} et/ou Y doivent être standardisées ou non, dans notre cas nous avons pris le partie de la version non standardisée, mais dans le cas contraire l’approche ne change pas.

La régression LASSO:

La régression LASSO (Least Absolute Shrinkage and Selection Operator), ou L_1-pénalisation, a été élaborée en 1995 par Robert Tibshirani. Elle est particulièrement adaptée au problème de type n <<<< P et notamment à la recherche d’un sous-ensemble de variables explicatives depuis le modèle complet.

La résolution du modèle LASSO est plus compliquée que celle du modèle classique et du modèle RIDGE puisqu’il est impossible d’obtenir les dérivés, solution du problème, étant donné qu’il s’agit d’un L_1-problème et donc non dérivable. Le problème peut être synthétisé de la manière suivante, estimer les paramètres \beta_L du modèle tel que,

\beta_L = argmin_{\beta \in R ^{P+1}} \lbrace ||Y - X \beta|| ^2 + \lambda \sum_{p = 1} ^P |\beta_p| \rbrace

, avec \lambda paramètre de L_1-pénalisation à déterminer.

add.png

Ce problème n’est résoluble que par un algorithme itératif. Plusieurs algorithmes de résolution du problème LASSO existent, le plus célèbre est l’algorithme LARS dont voici le déroulement.

Soit \lambda le paramètre de pénalisation à appliquer aux coefficients \beta du modèle construit par régression linéaire classique. s \in \lbrace -1, 0, +1 \rbrace ^P le vecteur des signes associés aux coefficients pénalisés \beta_{LARS}.

– Initialisation des paramètres: \lambda_0 = \infty, s_0 = c(0, \cdots, 0) et de taille 1 \times (P+1), \beta_{LARS} (\lambda_0) = \beta ^{s_0} (\lambda_0) = (0, \cdots, 0) et également de taille 1 \times (P+1).

– Première étape: diminuer \lambda_1 ^k \geq \lambda_0 jusqu’à ce que \exists p, k / |(X ^p) ^t (Y - \mathbf{X} \beta_{LARS})| > \lambda_1 ^k. Le paramètre \lambda_1 ^k correspondant à cette condition donne naissance à \lambda_1 = \lambda_1 ^k. Calculer ensuite son signe \epsilon_1 \in \lbrace -1 , +1 \rbrace, ce qui donne s_1 = s_0, s_1 ^p = \epsilon_1.

– Étapes suivantes (itération N° it \geq 2): diminuer \lambda_{it} ^k \geq \lambda_{it} jusqu’à ce que \exists p, k / |(X ^p) ^t (Y - \mathbf{X} \beta_{LARS})| > \lambda_{it} ^k, \beta_{LARS} = \beta ^{s_{it}} ou que \exists p, s_{it} ^p \ne 0. Définir \lambda_{it+1}, s^{it+1} à partir des résultats obtenus pour l’itération it.

– Fin de l’algorithme: dés que \lambda est nul c’est que nous avons enfin trouver le vecteur de coefficients \beta_{LARS} solution au problème LASSO.

Parmi les autres algorithmes les plus connus pour la résolution du problème LASSO, nous pouvons citer group-LASSO qui applique une pénalisation par sous-ensemble de variables et l’adaptative-LASSO qui agit par pondération de pénalisation directement sur \beta.

La régression LASSO implique que les variables explicatives avec un très faible pouvoir explicatif auront un coefficient nul, ce qui facilite la sélection des variables les plus intéressantes vis à vis de la variable réponse Y. Mais aussi le désavantage de centraliser l’information sur une seule des variables prise aléatoirement parmi toutes celles avec lesquelles elle est corrélée.

La régression ELASTIC NET:

La régression ELASTIC NET a été élaborée par Hui Zhou et Trevor Hastie en 2005 afin de pallier aux problèmes que peut rencontrer la régression linéaire standard et en s’appuyant sur les forces des régressions RIDGE et LASSO. Ainsi, la régression ELASTIC NET va tirer avantage de la capacité du premier à pouvoir attribuer aux coefficients des variables explicatives corrélées le même coefficient et de la capacité du second à pouvoir attribuer des coefficients nuls aux variables explicatives peu informatives.

La régression ELASTIC NET se base alors sur le modèle suivant,

\beta_E = argmin_{\beta} (||Y - \mathbf{X} \beta|| ^2 +  \lambda_{L_2} || \beta || ^2 + \lambda_{L_1} || \beta ||_1)

Nous retrouvons, respectivement, les normes L_1, L_2 des modèles respectifs LASSO et RIDGE ainsi que la partie commune inspirée de la régression linéaire.

add.png

De nombreux algorithmes ont été implémantés pour la résolution du problème ELASTIC NET, l’un des plus connus est celui de van der Kooij et nommé Coordinate Descent.

\bullet Annexe théorique:

– Nous présentons ici une esquisse de la démonstration de l’espérance et de la variance totale, entités qui sont à la base du modèle linéaire dont fait partie la régression linéaire.

Théorème de l’espérance totale: E[E[Y/X]] = E[Y]

En effet,

E[E[Y/X]] = \sum_x E[Y/X=x] P(X=x)

= \sum_x (\sum_y y P(Y=y/X=x)) P(X=x)

= \sum_y y \sum_x P(Y=y/X=x)P(X=x)

= \sum_y y P(Y=y)

= E[Y]

Théorème de la variance totale: V(Y) = E[V(Y/X)]+V(E[Y/X])

En effet,

V(Y) = E[Y - E[Y]] ^2

= E[Y - E[Y/X] + E[X/Y] - E[Y]]^2

= E[Y - E[Y/X]^2] + 2 E[(Y - E[Y/X]) (E[Y/X] - E[Y])] + E[E[Y/X] - E[Y]]^2

= E[V(Y/X)] + 0 + V(E[Y/X])

= E[V(Y/X)] + V(E[Y/X])

Le terme nul ce justifie par le fait que,

E[(Y - E[Y/X]) (E[Y/X] - E[Y])] = (E[Y/X] - E[Y])(E[Y - E[Y/X]]/X)

, avec E[Y/X] - E[Y] constant à X figé et (E[Y - E[Y/X]]/X) = 0 après développement.

– Nous présentons désormais les éléments clés qui permettent de justifier la méthode des moindres carrés pour l’estimations des coefficients \beta.

L’objectif que nous nous fixons est de trouver une équation de la forme,

\hat{Y} = \beta_0 1 + \beta_1 X ^1 + \cdots + \beta_P X ^P = \mathbf{X} \beta

, tel que ||Y - \hat{Y}|| ^2 soit minimal. Afin de pouvoir rendre cette contrainte vérifiable, nous partons du principe que nous nous situons dans R ^n et que nous nous munissons de la métrique \mathbf{D}. en outre, nous venons de définir le critère des moindres carrés.

\hat{Y} est alors la projection \mathbf{D}-orthogonale de Y sur le sous-espace W engendré par (1, X ^1, \dots, X ^P) et donc de dimension P + 1. Or, nous savons que l’opérateur \mathbf{D}-orthogonal sur W est de la forme,

\mathbf{X} (\mathbf{X} ^t \mathbf{D} \mathbf{X}) ^{-1} \mathbf{X} ^t \mathbf{D}

D’où,

\hat{Y} = \mathbf{X} (\mathbf{X} ^t \mathbf{D} \mathbf{X}) ^{-1} \mathbf{X} ^t \mathbf{D} Y

, et selon la définition ci-dessus de \hat{Y}, nous pouvons alors écrire,

\hat{Y} = \mathbf{X} \beta = \mathbf{X} (\mathbf{X} ^t \mathbf{D} \mathbf{X}) ^{-1} \mathbf{X} ^t \mathbf{D} Y

\Rightarrow \mathbf{X} \beta = \mathbf{X} (\mathbf{X} ^t \mathbf{D} \mathbf{X}) ^{-1} \mathbf{X} ^t \mathbf{D} Y

\Rightarrow \beta = (\mathbf{X} ^t \mathbf{D} \mathbf{X}) ^{-1} \mathbf{X} ^t \mathbf{D} Y

Enfin, si \mathbf{D} = \frac{1}{n} \mathbf{I}, alors,

\beta = (\mathbf{X} ^t \frac{\mathbf{I}}{n} \mathbf{X}) ^{-1} \mathbf{X} ^t \frac{\mathbf{I}}{n} Y = (\frac{\mathbf{I}}{n}) ^{-1} \frac{\mathbf{I}}{n} (\mathbf{X} ^t \mathbf{X}) ^{-1} \mathbf{X} ^t Y = (\mathbf{X} ^t \mathbf{X}) ^{-1} \mathbf{X} ^t Y

\bullet Exemple:

Soit le jeu de données ci-dessous,

add.png

Modélisons notre variable réponse Y à partir des variables explicatives X ^1, X ^2.

Régression linéaire:

Tout d’abord nous construisons \mathbf{X}_1 = (1, \mathbf{X}). Nous pouvons alors calculer les coefficients de le régression linéaire,

\beta = (\begin{pmatrix} 1 & 1 & \cdots & 1 & 1 \\ 8.1472 & 9.0579 & \cdots & 7.9221 & 9.5949 \\ 1.9593 & 2.5472 & \cdots & 11.5497 & 9.9172 \end{pmatrix} \times \begin{pmatrix} 1 & 8.1472 & 1.9593 \\ 1 & 9.0579 & 2.5472 \\ \cdots & \cdots & \cdots \\ 1 & 7.9221 & 11.5497 \\ 1 & 9.5949 & 9.9172 \\ \end{pmatrix}) ^{-1} \times \begin{pmatrix} 1 & 1 & \cdots & 1 & 1 \\ 8.1472 & 9.0579 & \cdots & 7.9221 & 9.5949 \\ 1.9593 & 2.5472 & \cdots & 11.5497 & 9.9172 \end{pmatrix} \times \begin{pmatrix} 0.8970 \\ 2.0949 \\ \cdots \\ 18.9751 & 19.8936 \\ \end{pmatrix}

= \begin{pmatrix} 0.43220908 & -0.03046141 & -0.0174284 \\ -0.03046141 & 0.004836861 & -0.00005.534307 \\ -0.01742184 & -0.00005534307 & 0.001659434 \\ \end{pmatrix} \times \begin{pmatrix} 209.7543 \\ 1395.1132 \\ 2640.6649 \\ \end{pmatrix}

= \begin{pmatrix} 2.1553479 \\ 0.2124134 \\ 0.6504938 \\ \end{pmatrix}

\Rightarrow \beta_0 = 2.1553479, \beta_1 = 0.2124134, \beta_2 = 0.6504938

\Rightarrow Y = 2.1553479 + 0.2124134 \times X ^1 + 0.6504938 \times X ^2

Calcul des indicateurs du modèle construit:

  • La vraisemblance,

Le vecteur des résidus, soit la différence entre prédictions et valeurs réelles de Y, est,

\epsilon = (-4.263435, -3.641405, -1.436032, \cdots, 7.166825, 7.623983, 9.249089)

\Rightarrow RMSE = \sqrt{\frac{\sum_{i = 1} ^{20} \epsilon_i ^2}{20-3}} = \sqrt{\frac{395.9338}{17}} = \sqrt{23.29022} = 4.825994

Nous pouvons dés lors calculer la vraisemblance du modèle qui nous permettra de calculer les différents indicateurs de performance.

L(\beta_0, \beta_1, \beta_2, Y_1, \cdots, Y_{20}) = \frac{1}{\sqrt{2 \pi 4.825994 ^2}}  \times e ^{- \frac{1}{2 \times 4.825994 ^2} \sum_{i = 1} ^{20} \epsilon_i  ^2}

= 0.1465204 \times e ^{- 0.02146824 \times 395.9338}

= 0.1465204 \times 0.000203468

= 0.00002981221

Nous devons également calculer la vraisemblance du modèle constant. Nous avons alors,

\epsilon_{Cst} = (-1.2583479, -0.0604479, 0.8753521, \cdots, 15.9976521, 16.8197521, 17.7382521)

\Rightarrow RMSE_{Cst} = \frac{2050.828}{17} = 120.6369

Et donc,

L(\beta_0, Y_1, \cdots, Y_{20}) = \frac{1}{\sqrt{2 \pi 120.6369 ^2}}  \times e ^{- \frac{1}{2 \times 120.6369 ^2} \sum_{i = 1} ^{20} (epsilon_{Cst})_i  ^2}

= 0.003306967 \times e ^{- 0.005861447 \times 2050.828}

= 0.003306967 \times 0.000006017615

= 0.00000001990005

  • Le R ^2,

Nous avons \overline{Y} = 10.48771, alors,

R ^2 = 1 - \frac{395.9338}{\sum_{i = 1} ^{20} e_i ^2 + \sum_{i = 1} ^{20} (\hat{Y}_i - 10.48771) ^2}

= 1 - \frac{395.9338}{395.9338 + 266.3272}

= 1 - \frac{395.9338}{662.261}

= 1 - 0.5978516

= 0.4021484

  • Le R_{adj} ^2,

Étant donné que notre modèle présente un coefficient constant, nous avons \sharp \lbrace \beta_0 \rbrace = 1. Alors,

R_{adj} ^2 = 1 - \frac{(20 - 1) \times (1 - 0.4021484)}{20 - 3} = 1 - \frac{19 \times 0.5978516}{17} = 1 - 0.6681871 = 0.3318129

  • L’AIC,

AIC = 20 \times ln(\frac{\sum_{i = 1} ^{20} \epsilon_i ^2}{20}) + 2 \times 3 = 20 \times 2.985515 + 6 = 65.7103

  • Le BIC,

Nous avons \sigma ^2 = var(\epsilon) = 20.83862, alors,

BIC = 20 \times ln(\frac{\sum_{i = 1} ^{20} \epsilon_i ^2}{20}) + 2 \times (2 + 3) \times \frac{20 \times 20.83862}{\sum_{i = 1} ^{20} \epsilon_i ^2} - 2 \times (\frac{20 \times 20.83862}{\sum_{i = 1} ^{20} \epsilon_i ^2}) ^2

= 20 \times 2.985515 + 10 \times 1.052631 - 2 \times 1.108033

= 59.7103 + 10.52631 - 2.216066

= 68.02054

  • Le C_p de Mallows,

C_p = \frac{\sum_{i = 1} ^{20} epsilon_i ^2}{20.83862} - 20 + 2 \times 2 = \frac{395.9338}{20.83862} - 20 + 4 = 19 - 16 = 3

Régression RIDGE:

Les coefficients estimés lors de la régression linéaire de base nous donne le vecteur,

\beta = (2.1553479, 0.2124134, 0.6504938)

Nous calculons l’écart-type de \epsilon_i = Y_i - \hat{Y_i}, \forall i \in [1,n],

\epsilon = (0.8970 - 2.1553479 - 0.2124134 \times X_1 ^1 - 0.6504938 \times X_1 ^2, \cdots, 19.8936 - 2.1553479 - 0.2124134 \times X_{20} ^1 - 0.6504938 \times X_{20} ^2)

= (0.8970 - 5.160435, 2.0949 - 5.736305, \cdots, 18.9751 - 11.351117, 19.8936 - 10.644511)

= (-4.263435, -3.641405, \cdots, 7.623983, 9.249089)

\Rightarrow sd_{\epsilon} = 4.564934

Par conséquent nous avons,

\lambda = \frac{2 \times 4.564934 ^2}{(\begin{pmatrix} 0.2124134 & 0.6504938 \\ \end{pmatrix}) \times \begin{pmatrix} 0.2124134 \\ 0.6504938 \\ \end{pmatrix}} = \frac{41.67724}{5.113786} = 8.149976

Nous pouvons désormais procéder à l’estimation des coefficients du modèle RIDGE,

\beta = (\begin{pmatrix} 1 & 1 & \cdots & 1 & 1 \\ 8.1472 & 9.0579 & \cdots & 7.9221 & 9.5949 \\ 1.9593 & 2.5472 & \cdots & 11.5497 & 9.9172 \end{pmatrix} \times \begin{pmatrix} 1 & 8.1472 & 1.9593 \\ 1 & 9.0579 & 2.5472 \\ \cdots & \cdots & \cdots \\ 1 & 7.9221 & 11.5497 \\ 1 & 9.5949 & 9.9172 \\ \end{pmatrix} - \begin{pmatrix} 8.149976 & 0 & 0 \\ 0 & 8.149976 & 0 \\ 0 & 0 & 8.149976 \\ \end{pmatrix}) ^{-1} \times \begin{pmatrix} 1 & 1 & \cdots & 1 & 1 \\ 8.1472 & 9.0579 & \cdots & 7.9221 & 9.5949 \\ 1.9593 & 2.5472 & \cdots & 11.5497 & 9.9172 \end{pmatrix} \times \begin{pmatrix} 0.8970 \\ 2.0949 \\ \cdots \\ 18.9751 & 19.8936 \\ \end{pmatrix}

= (\begin{pmatrix} 11.85002 & 128.4068 & 214.2557 \\ 128.40680 & 1023.0899 & 1382.4922 \\ 214.25570 & 1382.4922 & 2889.9703 \end{pmatrix}) ^{-1} \times \begin{pmatrix} 209.7543 \\ 1395.1132 \\ 2640.6649 \\ \end{pmatrix}

= \begin{pmatrix} -0.7697995 \\ 0.4197220 \\ 0.7700205 \\ \end{pmatrix}

\Rightarrow \beta_0 = -0.7697995, \beta_1 = 0.4197220, \beta_2 = 0.7700205

\Rightarrow Y = -0.7697995 + 0.4197220 \times X ^1 + 0.7700205 \times X ^2

\bullet Application informatique:

Procédure SAS: https://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#reg_toc.htm

Package et fonction R: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/lm.html

\bullet Bibliographie:

– Sur les degrés mesurés des méridiens, et sur les longueurs observées sur pendule de Pierre Simon de Laplace.

– Nouvelles méthodes pour la détermination des orbites des comètes de Adrien Marie Legendre.

– Statistique. Dictionnaire encyclopédique de Yadolah Dodge.

– The Elements of Statistical Learning de Trevor Hastie, Robert Tibshirani et Jerome Friedman.

– Data Mining et statistique décisionnelle. L’intelligence des données de Stéphane Tufféry.

– Probabilités, analyse des données et Statistique de Gilbert Saporta.

– Application of Rigde Analysis to Regression Problems de A. E. Hoerl.

– Regression shrinkage and selection via the lasso de Robert Tibshirani.

– Le Lasso, ou comment choisir parmi un grand nombre de variables à l’aide de peu d’observations d’Anisse Ismaili et de Pierre Gaillard.

– Le site: http://www.ats.ucla.edu/stat/mult_pkg/faq/general/Psuedo_RSquareds.htm

– Pratique de la Régression Linéaire Multiple. Diagnostic et sélection de variables de Ricco Rakotomalala.

– Regularization and Variable Selection via the Elastic Net de Hui Zhou et Trevor Hastie.

L’Analyse discriminante de Fisher

snedecor

\bullet Présentation:

L’Analyse discriminante de Fisher, qui a été élaborée en 1936 par Ronald Aylmer Fisher, permet de discriminer une variable réponse Y binaire ou polychotomique (K \geq 2 classes) à partir d’une matrice de P variables explicatives continues \mathbf{X} = (X ^1, \cdots, X ^P).

L’Analyse discriminante de Fisher est réputée pour être le classifieur optimale lorsque les données suivent une loi normale. Il permet également la construction de deux formes de frontière: linéaire ou quadratique, permettant ainsi de s’adapter à plusieurs type de distribution de la réponse Y. Ces deux versions de l’analyse discriminante de Fisher nécessite néanmoins des hypothèses d’utilisation solides, la normalité des données pour les deux cas plus celle d’homoscédasticité pour la cas linéaire.

Une version étendue de l’analyse discriminante de Fisher existe lorsque les données sont qualitatives: l’algorithme DISQUAL de Gilbert Saporta, qui consiste à faire une Analyse des Correspondances Multiples (ACM) et ensuite appliquer l’analyse discriminante de Fisher sur les axes factoriels.

Enfin, l’analyse discriminante de Fisher se base sur la distance de Mahalanobis-Fisher afin d’établir les coefficients de la règle décisionnelle.

\bullet L’analyse discriminante de Fisher:

Hypothèses préliminaires: Variables explicatives \mathbf{X} continues, variable réponse Y binaire ou polychotomique. Normalité et homoscédasticité des données pour le cas linéaire, normalité pour le cas quadratique.

L’Analyse discriminante de Fisher propose deux types de frontière, l’une de forme linéaire et l’autre de forme quadratique. L’idée de l’analyse discriminante de Fisher est de se baser sur l’étude des covariances intra et inter-classes ainsi que sur l’élaboration d’un classifieur bayésien optimum au sens de la probabilité de l’erreur. Cet outil est à mi-chemin entre la modélisation et la classification, impliquant qu’il faut en général passer par une méthode d’apprentissage pour valider la règle décisionnelle.

La fonction discriminante du classifieur bayésien:

Notons, \forall k \in [1,K], \mathbf{\Sigma}_k la matrice de variance-covariance et \mu_k le centre de gravité du nuage d’individus associé à la classe Cl_k de Y, soit à \mathbf{X|_{Y = Cl_k}}, d’effectif n_k.

Supposons que les observations de chaque classe soient générées selon une loi normale de paramètres (\mu_k, \mathbf{\Sigma}_k):

P(Y = Cl_k / \mathbf{X}) = \frac{1}{(2 \pi) ^{\frac{P}{2}} \vert \mathbf{\Sigma}_k \vert ^{\frac{1}{2}}} e^{-\frac{1}{2} (\mathbf{X} - \mu_k) ^T \mathbf{\Sigma}_k ^{-1} (\mathbf{X} - \mu_k)}

Si les matrices de variance-covariance sont identiques alors nous sommes dans le cas d’une discrimination linéaire (hypothèse de normalité et d’homoscédasticité), dans le cas inverse nous nous orientons vers du quadratique (hypothèse de normalité).

– Dans le premier cas, où \mathbf{\Sigma}_1 = \cdots = \mathbf{\Sigma}_K = \mathbf{\Sigma} (la matrice de variance-covariance toutes classes confondues), les fonctions discriminantes du classifieur bayésien deviennent:

\delta_k(\mathbf{X}) = (\mathbf{X} - \mu_k ^T) \cdot \mathbf{\Sigma} ^{-1} \cdot (\mathbf{X} - \mu_k) - 2 log(\frac{n_k}{n})

– Dans le second cas, où \mathbf{\Sigma}_1 \neq \cdots \neq \mathbf{\Sigma}_K, elles valent:

\delta_k(\mathbf{X}) = (\mathbf{X} - \mu_k ^T) \cdot \mathbf{\Sigma}_k ^{-1} \cdot (\mathbf{X} - \mu_k) + log(\mid \mathbf{\Sigma}_k \mid) - 2 log(\frac{n_k}{n})

Le calcul des coefficients:

Le vecteur des coefficients associé à la règle décisionnelle se retrouve au travers de la formule de la fonction discriminante décrite ci-dessus en fonction des deux versions,

– pour la règle linéaire, \beta = (\mu_k ^T \mathbf{\Sigma}_k ^{-1} \mu_k, \mathbf{\Sigma}_k \mu_k)

– pour la règle quadratique, \beta = (\mu_k ^T \mathbf{\Sigma} ^{-1} \mu_k, \mathbf{\Sigma} \mu_k)

Le premier terme du vecteur représentant le coefficient constant.

Règle décisionnelle:

Nous cherchons à savoir de quelle classe Y = Cl_k le nouvel individu i que nous voulons classer est le plus près. Nous répondons à cette interrogation via la distance définie par la matrice \mathbf{\Sigma}_k et qui correspond la règle descriptive de Mahalanobis-Fisher, de définition,

d^2 (X_i,\mu_k) = X_i ^T \cdot \mathbf{\Sigma}_k ^{-1} \cdot X_i - 2 \cdot (X_i \cdot (\beta_2, \cdots, \beta_K) - \beta_1) + C

Où,

C = 0 dans le cas linéaire,

C = log( | \mathbf{\Sigma}_k |) - 2 \cdot log(\frac{n_k}{n}) dans le cas quadratique.

La règle décisionnelle basée sur le théorème de Bayes nous donne:

P(Y = Cl_k / X_i) = \frac{P_k \cdot P(X_i / Cl_k)}{\sum_{k = 1} ^K P_k \cdot P(X / Cl_k)} = \frac{P_k \cdot e^{-\frac{1}{2} \cdot d ^2(X_i, \mu_k)}}{\sum_{k = 1} ^K P_k \cdot e^{-\frac{1}{2} \cdot d ^2(X_i, \mu_k)}}

, avec P_k = \frac{n_k}{n}.

\bullet Annexe théorique:

Nous présentons ici une esquisse de la démonstration de la fonction de Mahalanobis-Fisher.

Pour P variables explicatives, nous pouvons rechercher la combinaison linéaire définie par les coefficients \beta qui donnent la valeur maximale pour la statistique de test,

\frac{\beta ^T \mathbf{B} \beta}{\beta ^T \mathbf{\Sigma} \beta}

, où \mathbf{B} est la matrice de variance des K centres de gravité.

La solution est donnée par,

\mathbf{W} ^{-1} \mathbf{B} \beta = \delta \beta

, avec \delta maximal et \mathbf{W} la matrice de variance intra-classe.

Les vecteurs propres de \mathbf{W} ^{-1} \mathbf{B} sont les mêmes que ceux de \mathbf{\Sigma} ^{-1} \mathbf{B} avec \delta = \frac{\lambda}{1 - \lambda}.

En effet,

\mathbf{B} \beta = \lambda (\mathbf{W} + \mathbf{B}) \delta \Rightarrow (1 - \lambda) \mathbf{B} \beta = \lambda \mathbf{W} \beta

D’où,

\mathbf{W} ^{-1} \mathbf{B} \beta = \frac{\lambda \beta}{1 - \lambda}

L’utilisation de \mathbf{\Sigma} ^{-1}, \mathbf{W} ^{-1} est donc indifférent. La seconde porte le nom de métrique de Mahalanobis.

\bullet Exemple:

Soit la base de données suivante,

add

Avec « Statut » notre variable réponse et (X ^1, X ^2) nos variables explicatives.

Notons que nous avons n_1 = n_2 = 10 \Rightarrow P_1 = P_2 = \frac{10}{20} = \frac{1}{2}.

Nous calculons pour le groupe des \mathbf{X}|_{Y = 1} la matrice de variance-covariance,

\mathbf{\Sigma}_1 = \begin{pmatrix} 11.961355 & 1.749008 \\ 1.749008 & 2.129509 \end{pmatrix}

, de déterminant \mid \mathbf{\Sigma}_1 \mid = 22.41278. Et le barycentre,

\mu_1 = (6.23856, 4.99603)

Et pour le groupe des \mathbf{X}|_{Y = 2},

\mathbf{\Sigma}_2 = \begin{pmatrix} 10.9457223 & 0.3020457 \\ 0.3020457 & 15.7377368 \end{pmatrix}

, de déterminant \mid \mathbf{\Sigma_2} \mid = 172.1697. Et le barycentre,

\mu_2 = (6.60212, 5.52999)

Pour le cas quadratique:

Nous obtenons alors les vecteurs de coefficients:

\beta_1 = \frac{1}{2} \cdot (\mu_1 ^T \cdot \mathbf{\Sigma}_1 ^{-1} \cdot \mu_1, \mathbf{\Sigma}_1 ^{-1} \cdot \mu_1) = (6.07717, 0.2028741, 2.1794703)

\beta_2 = \frac{1}{2} \cdot (\mu_2 ^T \cdot \mathbf{\Sigma}_2 ^{-1} \cdot \mu_2, \mathbf{\Sigma}_2 ^{-1} \cdot \mu_2) = (2.900191, 0.5937870, 0.3399878)

En se reportant à la formule de la distance de Mahalanobis-Fisher nous déterminons les quantités d_1 ^2 et d_2 ^2 pour l’individu N°1 et de caractéristiques e = (8.1472, 3.1101):

d_1 ^2 = e ^T \times \mathbf{\Sigma}_1 ^{-1} \times e - 2 \times [ e \times (0.2028741, 2.1794703) ^T - 6.07717] + log(22.41278) -2 \times log(\frac{10}{20})

 = 7.5142 - 4.708112 + 3.109631 - (-1.386294)

= 7.302015

d_2 ^2 = e ^T \times \mathbf{\Sigma}_2 ^{-1} \times e - 2 \times [ e \times (0.5937870, 0.3399878) ^T - 2.900191 ] + log(172.1697) -2 \times log(\frac{10}{20})

 = 6.593437 - 5.989814 + 5.148481 - (-1.386294)

= 7.138398

Nous pouvons désormais calculer sa probabilité d’appartenir à la classe Y = 1:

P (Y =1 \setminus (8.1472, 3.1101)) = \frac{e ^{- \frac{1}{2} \cdot 7.302015}}{e ^{- \frac{1}{2} \cdot 7.302015} + e ^{\frac{1}{2} \cdot 7.138398}} = 0.4795593

, soit 48 \% de chance que cet individu soit Y = 1, donc nous le classerons Y = 2.

En effectuant ce calcul pour tous les individus nous classons chacun des vingt individus. Nous obtenons la matrice de confusion suivante:

\begin{tabular}{|l|c|c|} \hline & Predit Y = 1 & Predit Y = 2 \\ \hline Y = 1 & 9 & 1 \\ \hline Y = 2 & 0 & 10 \\ \hline \end{tabular}

Pour le cas linéaire:

Nous devons calculer la matrice de covariance commune à toute la population:

\mathbf{\Sigma} = \frac{20}{19} \times \begin{pmatrix} 10.885504 & 1.022638 \\ 1.022638 & 8.538462 \end{pmatrix} = \begin{pmatrix} 11.458425 & 1.076461 \\ 1.076461 & 8.987855 \end{pmatrix}

Analoguement à la construction de la forme quadratique, nous avons les vecteurs des coefficients qui valent:

\beta_1 = \frac{1}{2} \cdot (\mu_1 ^T \cdot \mathbf{\Sigma} ^{-1} \cdot \mu_1, \mathbf{\Sigma} ^{-1} \cdot \mu_1) = (2.7924940, 0.4978326, 0.4962400)

\beta_2 = \frac{1}{2} \cdot (\mu_2 ^T \cdot \mathbf{\Sigma} ^{-1} \cdot \mu_2, \mathbf{\Sigma} ^{-1} \cdot \mu_2) = (3.2582809, 0.5242776, 0.5524818)

Enfin, les quantités d_1 ^2 et d_2 ^2 pour le même individu (de caractéristiques: e = (8.1472, 3.1101)):

d_1 ^2 = e ^T \times S_{Y_1} ^{-1} \times e - 2 \times [ e \times (0.4978326, 0.4962400) ^T - 2.7924940 ]

= 6.411483 - 5.613607

= 0.797876

d_2 ^2 = e ^T \times S_{Y_2} ^{-1} \times e - 2 \times [ e \times (0.5242776, 0.5524818) ^T - 3.2582809 ]

= 6.411483 - 5.462773

= 0.94871

Nous pouvons désormais calculer sa probabilité d’appartenir à la classe Y_1:

P (Y = 1 \setminus (8.1472, 3.1101)) = \frac{e ^{- \frac{1}{2} \cdot 0.797876}}{e ^{- \frac{1}{2} \cdot 0.797876} + e ^{\frac{1}{2} \cdot 0.94871}} = 0.5188452

, soit 52\% de chance que cet individu soit Y_1, donc nous le classerons Y_1.

\begin{tabular}{|l|c|c|} \hline & Predit Y1 & Predit Y2 \\ \hline Y1 & 5 & 5 \\ \hline Y2 & 5 & 5 \\ \hline \end{tabular}

\bullet Application informatique:

Procédure SAS: https://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_discrim_sect004.htm

Package et fonction R

https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/lda.html

https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/qda.html

\bullet Bibliographie:

– Probabilité, analyse des données et Statistique de Gilbert Saporta

– Data mining et statistique décisionnelle. L’intelligence des données de Stéphane tufféry

– Comprendre et utiliser les statistiques dans les sciences de la vie de Bruno Falissard

– The elements of statistical learning de Trevor Hastie, Robert Tibshirani et Jerome Friedman

Les K-plus proches voisins pondérés

add2.png

\bullet Présentation:

Les K-plus proches voisins pondérés, publié en 1951 par Evelyn Fix et Joseph L. Hodges et noté KKNN, est un outil permettant de discriminer une variable réponse Y binaire ou polychotomique (C \geq 2 classes) à partir d’une matrice de P variables explicatives continues \mathbf{X} = (X ^1, \cdots, X ^p).

L’Algorithme des K-plus proches voisins pondérés fait partie des outils de discrimination par apprentissage statistique. Il permet notamment d’élaborer des frontières de toutes les formes possibles et offre l’avantage de ne pas nécessiter d’hypothèse à priori sur la distribution de \mathbf{X}.

L’Idée de l’algorithme est de déterminer la classe d’une observation en fonction d’une distance pondérée, par une fonction noyau, entre elle et ses plus proches voisins. Par conséquent, les K plus proches voisins pondérés nécessitent de fixer trois paramètres afin de pouvoir être utilisé.

\bullet Les K plus proches voisins pondérés:

Hypothèse préliminaire: Variables explicatives \mathbf{X} continues, variable réponse Y binaire ou polychotomique.

L’Algorithme des K-plus proches voisins pondérées nécessite de fixer trois paramètres: la fonction noyau f (pondération des distances, que nous pouvons retrouver sous l’appellation kernel dans la littérature), la fonction distance d et le nombre de voisins K. Le choix des paramètres optimaux se fait selon le critère de maximisation des performances.

La fonction noyau f:

Le vote des plus proches voisins est équivalent au mode de distribution de la classe. L’idée étant de pondérer la distance entre l’observation ciblée et celles l’entourant.

Les huit fonctions noyaux suivantes sont les plus utilisées pour l’algorithme des K-plus proches voisins pondérés.

– Rectangulaire: f(d) = \frac{1}{2} \cdot I(\vert d \vert \leq 1)

– Triangulaire: f(d) = (1 - \vert d \vert) \cdot I(\vert d \vert \leq 1)

– Epanechinikov: f(d) = \frac{3}{4} \cdot (1 - d^2) \cdot I(\vert d \vert \leq 1)

– Bi-poids: f(d) = \frac{15}{16} \cdot (1 - d^2) ^2 \cdot I(\vert d \vert \leq 1)

– Tri-poids: f(d) = \frac{35}{32} \cdot (1 - d^2) ^3 \cdot I(\vert d \vert \leq 1)

– Cosinus: f(d) = \frac{\pi}{4} \cdot cos(\frac{\pi}{2} d) \cdot I(\vert d \vert \leq 1)

– Gaussien: f(d) = \frac{1}{\sqrt{2 \pi}} \cdot e ^{-\frac{1}{2} d^2}

– Inverse: f(d) = \frac{1}{\vert d \vert}

Le graphe ci-dessous présente la densité des différents noyaux présentés:

add

Le choix de la distance d:

Les six distances principales utilisées afin de calculer la matrice de distance entre les différents voisins sont:

– La distance euclidienne: d(X_{i_1}, X_{i_2}) = \sqrt{\sum_{j = 1} ^P (X_{i_1} ^j - X_{i_2} ^j) ^2}

– La distance maximale: d(X_{i_1}, X_{i_2}) = max_{j \in [1, P]} | X_{i_1} ^j - X_{i_2} ^j |

– La distance de Manhattan: d(X_{i_1}, X_{i_2}) = \sum_{j = 1} ^P | X_{i_1} ^j - X_{i_2} ^j |

– La distance de Canberra: d(X_{i_1}, X_{i_2}) = \sum_{j = 1} ^P \frac{| X_{i_1} ^j - X_{i_2} ^j |}{| X_{i_1} ^j | + | X_{i_2} ^j |}

– La distance de Minkowski (pour q, paramètre à fixer), qui est entre autre une généralisation de la distance euclidienne: d(X_{i_1}, X_{i_2}) = (\sum_{j = 1} ^P (X_{i_1} ^j - X_{i_2} ^j) ^q) ^{\frac{1}{q}}

L’algorithme KKNN:

Concluons la présentation de la méthodologie des K-plus proches voisins pondérés par l’algorithme d’élaboration de la règle décisionnelle.

– Étape 1: définition du nombre de voisins K, de la fonction distance d et de la fonction noyau f.

– Étape 2: Soit l’échantillon d’apprentissage L de définition  \lbrace (Y_i, X_i), i \in [1,n_L]. Pour l’ensemble des observations i \in [1,n_L] nous appliquons les étapes 3 à 6 suivantes.

– Étape 3: Sélection des (K + 1) plus proche voisins de X_i selon d(X,X_i), \forall i \in [1,n_L].

– Étape 4: Standardisation des K plus proches petites distances via le (K + 1)^{eme} plus proche voisin afin d’obtenir des valeurs comprises entre [0,1]:

D(i) = D(X,X_i) = \frac{d(X,X_i)}{d(X,X_k)} \forall i \in [1,K]

– Étape 5: Transformation des distances normalisées D(i) en pondération w(i) à partir de la fonction noyau:

w(i) = f(D(i))

– Étape 6: La classe de X_i est déterminée via vote majoritaire au sein des K plus proches voisins pour les différentes classes C de Y:

\overbrace{y} = max_{c \ in [1,C]} (\frac{1}{\sharp \lbrace \mbox {Voisins consideres de classe Y = c} \rbrace} \times \sum_{i = 1} ^k w(i) \cdot I(y(i) = c))

– Étape 7: Une fois les étapes 3 à 6 parcourues pour l’ensemble des observations i \in [1,n_L], nous obtenons une carte des secteurs des différentes classes de Y à partir de laquelle nous pouvons prédire la classe d’une nouvelle observation. Le modèle est alors vérifier avec le jeu test.

\bullet Annexe théorique:

Nous rappelons dans cette section de l’article les propriétés qu’une fonction f doit rempir pour être une fonction noyau,

\int_{- \infty} ^{+ \infty} f(x) dx = 1

f(d) \geq 0, \forall d \in R

f(d) \rightarrow max pour d = 0

f(d) \searrow quand d \rightarrow + \infty

\bullet Exemple:

Soit la base de données suivante,

add

Avec « Statut » notre variable réponse et (X ^1, X ^2) nos variables explicatives. Fixons les paramètres suivants pour l’étape 1:

– Nombre de voisins à considérer: K = 3,

– Le kernel sera la fonction noyau inverse,

– La fonction distance d est la distance de Minkowski pour q = 1.

L’étape 2, qui consiste au calcul de la matrice des distances via le paramètre d, nous donne:

\begin{pmatrix} & i_1 & i_2 & i_3 & i_4 & i_5 & i_6 & i_7 & i_8 & \ldots \\ i_2 & 1.9014 & . & . & . & . & . & . & . & \ldots \\ i_3 & 8.5548 & 8.4748 & . & . & . & . & . & . & \ldots \\ i_4 & 4.9442 & 3.0428 & 10.1440 & . & . & . & . & . & \ldots \\ i_5 & 4.7993 & 4.7193 & 6.3519 & 3.7921 & . & . & . & . & \ldots \\ i_6 & 8.9926 & 8.9126 & 0.4378 & 10.2952 & 6.5031 & . & . & . & \ldots \\ i_7 & 6.2970 & 6.3288 & 2.2578 & 9.3716 & 5.5795 & 2.6956 & . & . & \ldots \\ i_8 & 2.7784 & 4.6798 & 5.9764 & 7.7226 & 3.9305 & 6.4142 & 3.7186 & . & \ldots \\ i_9 & 4.2673 & 2.3659 & 9.4671 & 1.5595 & 3.3878 & 9.6183 & 8.6947 & 7.0457 & \ldots \\ i_{10} & 5.2645 & 3.3631 & 10.4643 & 0.7099 & 4.1124 & 10.6155 & 9.6919 & 8.0429 & \ldots \\ i_{11} & 8.5914 & 10.4928 & 4.0040 & 13.5356 & 9.7435 & 4.4418 & 4.1640 & 5.8130 & \ldots \\ i_{12} & 2.6820 & 2.7620 & 11.2368 & 5.6530 & 7.4813 & 11.6746 & 8.9790 & 5.2604 & \ldots \\ i_{13} & 1.5493 & 1.6293 & 10.1041 & 4.5203 & 6.3486 & 10.5419 & 7.8463 & 4.1277 & \ldots \\ i_{14} & 10.1913 & 10.1113 & 8.8043 & 7.2203 & 5.3920 & 8.9555 & 8.0319 & 7.6129 & \ldots \\ i_{15} & 5.9395 & 5.8595 & 10.8505 & 2.9685 & 4.4986 & 11.0017 & 10.0781 & 8.4291 & \ldots \\ i_{16} & 11.6593 & 11.5793 & 3.4025 & 8.6883 & 6.8600 & 3.5537 & 5.3623 & 9.0809 & \ldots \\ i_{17} & 4.9571 & 6.8585 & 5.6527 & 9.9013 & 6.1092 & 6.0905 & 3.3949 & 2.1787 & \ldots \\ i_{18} & 3.0667 & 3.1467 & 11.6215 & 6.0377 & 7.8660 & 12.0593 & 9.3637 & 5.6451 & \ldots \\ i_{19} & 6.1799 & 6.0999 & 10.9295 & 3.2089 & 4.5776 & 11.0807 & 10.1571 & 8.5081 & \ldots \\ i_{20} & 8.4202 & 6.5188 & 13.6200 & 3.4760 & 7.2681 & 13.7712 & 12.8476 & 11.1986 & \ldots \\ \end{pmatrix}

Nous cherchons à classer le premier individu i_1. Nous commençons donc par la recherche des trois voisins les plus proches. En regardant dans la table ci-dessus, nous trouvons qu’il s’agit des observations i_{13}, i_2 et i_{12}. Enfin, le (K+1) ^{ieme} = 4 ^{ieme} plus proche voisin qui servira à la standardisation des distances est i_8.

En étape 3, nous standardisons donc les distances entre i_1 et ces trois observations relevés par rapport à sa distance avec i_8:

D_{i_1} (i_{13}) = \frac{d(X_{i_1}, X_{i_{13}})}{d(X_{i_1}, X_{i_8})} = \frac{1.5493}{2.7784} = 0.5576231

D_{i_1} (i_2) = \frac{d(X_{i_1}, X_{i_2})}{d(X_{i_1}, X_{i_8})} = \frac{1.9014}{2.7784} = 0.6843507

D_{i_1} (i_{12}) = \frac{d(X_{i_1}, X_{i_{12}})}{d(X_{i_1}, X_{i_8})} = \frac{2.6820}{2.7784} = 0.9653038

L’avant dernière étape pour la classification de l’observation i_1 nous amène à la pondération en fonction de la fonction noyau f paramétrée. Nous obtenons alors:

w_{i_1} ( i_{13} ) = \frac{1}{\vert 0.5576231 \vert} = 1.793326

w_{i_1} ( i_2 ) = \frac{1}{\vert 0.6843507 \vert} = 1.461239

w_{i_1} ( i_{12} ) = \frac{1}{\vert 0.9653038 \vert} = 1.035943

Pour l’étape 5, nous allons évaluer les deux classifications possibles et voir pour laquelle nous obtenons le plus de voix.

Nous avons,

V_{Y = 1} = \frac{1}{\sharp \lbrace \mbox{ Voisins consideres de classe Y = 1 } \rbrace} \times \sum_{i = 1} ^k w_{i_1}(i) \cdot I(y(i) = 1)

 = \frac{1}{1} \times (1.793326 \times 0 + 1.461239 \times 1 + 1.035943 \times 0)

 = 1.461239

Ensuite, nous avons:

V_{Y = 2} = \frac{1}{\sharp \lbrace \mbox{ Voisins consideres de classe Y = 2 } \rbrace} \times \sum_{i = 1} ^k w_{i_1}(i) \cdot I(y(i) = 2)

 = \frac{1}{2} \times (1.793326 \times 1 + 1.461239 \times 0 + 1.035943 \times 1)

 = 1.414635

Nous constatons que V_{Y = 2} = 1.414635 < 1.461239 = V_{Y = 1} soit que nous classons i_1 chez les patients malades.

En appliquant ce raisonnement à toutes les observations nous obtenons la matrice de confusion et la carte des individus suivantes:

\begin{tabular}{|l|c|c|} \hline & Predit Y1 & Predit Y2 \\ \hline Y1 & 10 & 0 \\ \hline Y2 & 4 & 6 \\ \hline \end{tabular}

\bullet Application informatique:

Procédure SAS:

%macro WKNN(DATA =, Y =, d = , k = , kernel =);

/* La macro suivante permet d’appliquer l’algorithme des K-plus proches voisins pondérés pour P variables et Y binaire ou polychotomique
Les paramètres nécessaires sont: – DATA la base de données
– Y la variable réponse
– d la distance à appliquer et dont le codage peut être retrouver sur le site de la proc DISTANCE de SAS
– k le nombre de voisins
– kernel le noyau (Rectangulaire, Triangulaire, Epanechinikov, Bipoids, Tripoids, Cosinus, Gaussien, Inverse)
Le programme renvoie la base de données DATA avec ses prédictions selon l’outil WKNN ainsi que la matrice de confusion */

/* Initialisation des options pour épurer la log */
options nonotes spool;

/* Récupération du nombre d’observations */
proc sql noprint;
select count(*) into: n from &DATA.;
quit;

/* Sortie de la liste des variables explicatives */
proc contents data = &DATA. (drop = &Y.) out = biblio noprint;
run;

proc sql noprint;
select distinct(name) into: list_vars separated by  »  » from biblio;
quit;

/* Calcul de la matrice de distance selon la distance d paramétrée */
proc distance data = &DATA. method = &d. out = distX;
var interval(&list_vars.);
run;

/* Création de l’id unique */
data distX;
set distX;
id = _n_;
run;

/* Pour la récupération du statut réel, du nombre de classes et de leur nom respectif */
data Y;
set &DATA.;
Obs = _n_;
keep &Y. Obs;
run;

proc freq data = Y noprint;
tables &Y. / out = list_class;
run;

proc sql noprint;
select distinct(&Y.) into: list_class separated by  »  » from list_class;
select count(*) into: nb_class from list_class;
quit;

%let nmoinsun = %eval(&n.-1);

/* La matrice de distance produite par la proc distance n’est pas symétrique, nous menons les manipulations pour quel le soit */
%do i1 = 1 %to &nmoinsun.;

%let it = %eval(&i1. + 1);

%do i2 = &it. %to &n.;

data _null_;
set distX;
keep Dist&i1.;
if id = &i2.;
call symput(« valeur »,Dist&i1.);
run;

data distX;
set distX;
if id = &i1. then Dist&i2. = &valeur.;
run;

%end;

%end;

%let voisins = %eval(&k. + 1);

/* Détermination pour chaque observation de la classe selon l’outil WKNN */
%do i = 1 %to &n.;

data plusprochesvoisins;
run;

/* Récupération des distances entre l’observation i et les autres */
data distXencours;
set distX;
keep Dist&i. id;
if id ne &i.;
run;

/* Recherche des k plus proches voisins selon d */
%do v = 1 %to &voisins.;

proc sql noprint;
create table min as select min(Dist&i.) as min, Dist&i., id from DistXencours;
quit;

data min;
set min;
if Dist&i. = min;
run;

data min;
set min;
if _n_ = 1;
call symput(« min »,min);
call symput(« obs »,id);
run;

/* Nous complétons la table de synthèse des k plus proches voisins ainsi que du k+1 afin de pouvoir standardiser les calculs */
data plusprochesvoisins;
set plusprochesvoisins end = eof;
output;
if eof then do;
Obs = &obs.;
Dist = &min.;
ordre = &v.;
output;
end;
run;

data DistXencours;
set DistXencours;
if id ne &obs.;
run;

%end;

/* Calcul des pondérations selon le kernel paramétré */
data _null_;
set plusprochesvoisins;
if ordre = &voisins.;
call symput(« Denom »,Dist);
run;

data plusprochesvoisins;
set plusprochesvoisins;
if Obs ne .;
Std = Dist/&Denom.;
%if &kernel. = Rectangulaire %then %do;
if Dist < 1 or Dist = 1 then W = 0.5; if Dist > 1 then W = 0;
%end;
%if &kernel. = Triangulaire %then %do;
if Dist < 1 or Dist = 1 then W = 1 – abs(Dist); if Dist > 1 then W = 0;
%end;
%if &kernel. = Epanechinikov %then %do;
if Dist < 1 or Dist = 1 then W = 0.75 * (1 – Dist **2); if Dist > 1 then W = 0;
%end;
%if &kernel. = Bipoids %then %do;
if Dist < 1 or Dist = 1 then W = (15/16) * (1 – Dist **2) **2; if Dist > 1 then W = 0;
%end;
%if &kernel. = Tripoids %then %do;
if Dist < 1 or Dist = 1 then W = (35/32) * (1 – Dist **2) **3; if Dist > 1 then W = 0;
%end;
%if &kernel. = Cosinus %then %do;
if Dist < 1 or Dist = 1 then W = (constant(« pi »)/4) * cos(Dist*constant(« pi »)/2); if Dist > 1 then W = 0;
%end;
%if &kernel. = Gaussien %then %do;
W = (1/sqrt(2*constant(« pi »))) * exp(-0.5 * Dist **2);
%end;
%if &kernel. = Inverse %then %do;
W = 1/abs(Std);
%end;
run;

proc sort data = plusprochesvoisins;
by Obs;
run;

/* Fusion entre les calculs fait et le statut réel pour le vote majoritaire */
data plusprochesvoisins;
merge plusprochesvoisins Y;
by Obs;
if Dist ne .;
run;

proc freq data = plusprochesvoisins noprint;
tables &Y. / nopercent out = eff_Y;
run;

/* Calcul du score pour chaque classe en fonction des k voisins retenus */
%do c = 1 %to &nb_class.;

data eff_Yc;
set eff_Y;
if &Y. = « %scan(&list_class.,&c.) »;
call symput(« eff_c »,count);
run;

proc sql noprint;
select count(*) into: verif from eff_Yc;
quit;

%if &verif. = 0 %then %do;
%let eff_c = 0;
%end;

data calc_vote;
set plusprochesvoisins;
if &Y. ne « %scan(&list_class.,&c.) » then W = 0;
run;

proc sql noprint;
select sum(W) into: somme from calc_vote;
quit;

data plusprochesvoisins;
set plusprochesvoisins;
%if &eff_c. > 0 %then %do;
V_&c. = (1/&eff_c.) * &somme.;
%end;
%if &eff_c. = 0 %then %do;
V_&c. = 0;
%end;
run;

%end;

/* Choix de la classe de prédiction en fonction du vote */
proc contents data = plusprochesvoisins out = list_V noprint;
run;

data list_V;
set list_V;
if index(name, »V_ ») > 0;
run;

proc sql noprint;
select distinct(name) into: list_vote separated by « , » from list_V;
select distinct(name) into: list_vote2 separated by  »  » from list_V;
quit;

data plusprochesvoisins;
set plusprochesvoisins;
maxV = max(&list_vote.);
call symput(« max »,maxV);
run;

proc transpose data = plusprochesvoisins out = prediction;
var &list_vote2.;
run;

data prediction;
set prediction;
max = &max.;
compar_i = put(COL1,10.);
compar_v = put(max,10.);
if compar_i = compar_v then pred = _n_;
keep _name_ compar_i compar_v pred;
call symput(« pred »,pred);
run;

data _null_;
set prediction;
if pred ne .;
call symput(« pred »,pred);
run;

data &DATA.;
set &DATA.;
if _n_ = &i. then Pred_WKNN = &Pred.;
run;

%put Observation &i. traitée;

%end;

/* Sortie de la matrice de confusion dans la log*/
proc freq data = &DATA.;
tables &Y. * Pred_WKNN;
run;

/* Suppression des tables temporaires et inutiles */
proc datasets lib = work nolist;
delete distX biblio Y list_class plusprochesvoisins distXencours min eff_Y eff_Yc calc_vote list_V prediction;
run;

/* Reset des options d’affichage de la log */
options notes nospool;

%mend;

Package et fonction R: http://www.inside-r.org/packages/cran/kknn/docs/kknn

\bullet Bibliographie:

– An important contribution to nonparametric discriminant analysis and density estimation de Evelyn Fix et Joseph L. Hodges

– Algorithme des K-plus proches voisins pondérés et application en diagnostic de Eve Mathieu-Dupas

– The elements of statistical learning de Trevor Hastie, Robert Tibshirani et Jerome Friedman

– Theory of reproducing Kernel de N. Aronszajn