Econométrie et Machine Learning
Arthur Charpentier Emmanuel Flachaire Antoine Ly
stat.OT econ.EM














Économétrie & Machine Learning

















Arthur Charpentier

Université de Rennes 1 & CREM

7 Place Hoche, 35065 Rennes Cedex, France

arthur.charpentier@univ-rennes1.fr








Emmanuel Flachaire

Aix-Marseille Université, AMSE, CNRS & EHESS

5 bd Maurice Bourdet, CS 50498, 13205 Marseille Cedex 01, France

emmanuel.flachaire@univ-amu.fr




















et




Antoine Ly

Université Paris-Est

5, boulevard Descartes, 77454 Marne-la-Vallée cedex, France

antoine.ly@gmail.com 



















L'économétrie et l'apprentissage machine semblent avoir une finalité en commun: construire un modèle prédictif, pour une variable d'intérêt, à l'aide de variables explicatives (ou features). Pourtant, ces deux champs se sont développés en parallèle, créant ainsi deux cultures différentes, pour paraphraser Breiman. Le premier visait à construire des modèles probabilistes permettant de décrire des phénomèmes économiques. Le second utilise des algorithmes qui vont apprendre de leurs erreurs, dans le but, le plus souvent de classer (des sons, des images, etc). Or récemment, les modèles d'apprentissage se sont montrés plus efficaces que les techniques économétriques traditionnelles (avec comme prix à payer un moindre pouvoir explicatif), et surtout, ils arrivent à gérer des données beaucoup plus volumineuses. Dans ce contexte, il devient nécessaire que les économètres comprennent ce que sont ces deux cultures, ce qui les oppose et surtout ce qui les rapproche, afin de s'approprier des outils développés par la communauté de l'apprentissage statistique, pour les intégrer dans des modèles économétriques.



JEL Code: C18; C52; C55

Key-words: apprentissage; données massives; économétrie; modélisation; moindres carrés;





Mars  2018














§ INTRODUCTION


L'utilisation de techniques quantitatives en économie remonte probablement au 16ème siècle, comme le montre Morgan. Mais il faudra attendre le début du XXième siècle pour que le terme « économétrie » soit utilisé pour la première fois, donnant naissance à l'Econometric Society en 1933. Les techniques de machine learning (apprentissage machine) sont plus récentes. C'est à Arthur Samuel, considéré comme le père du premier programme d'auto-apprentissage, que l'on doit le terme « machine learning » qu'il définit comme « a field of study that gives computer the ability without being explicitly programmed ». Parmi les premières techniques, on peut penser à la théorie des assemblées de neurones proposée dans Hebb (qui donnera naissance au perceptron dans les années 1950, puis aux réseaux de neurones) dont Widrow montreront quinze ans plus tard les liens avec les méthodes des moindres carrés, aux SVM (support vector machine) et plus récemment aux méthodes de boosting. Si les deux communautés ont grandi en parallèle, les données massives imposent de créer des passerelles entre les deux approches, en rapprochant les « deux cultures » évoquées par Breiman, opposant la statistique mathématique (que l'on peut rapprocher de l'économétrie traditionnelle, comme le note Aldrich) à la statistique computationnelle, et à l'apprentissage machine de manière générale.



 §.§ La Modélisation économétrique


L'économétrie et les techniques d'apprentissage statistique supervisé sont proches, tout en étant très différentes. Proches au départ, car toutes les deux utilisent une base (ou un tableau) de données, c'est à dire des observations { (y_i,_i) }, avec i=1,⋯,n, _i∈𝒳⊂ℝ^p et y_i∈𝒴. Si y_i est qualitative, on parlera d'un problème de classification[Nous utiliserons ici le terme « classification » lorsque 𝒴 est un ensemble de classes, typiquement une classification binaire, 𝒴={0,1}, ce cas correspondant à la réalisation d'une variable indicatrice, 1_Y_t≤ 0, ou 1_Y∈𝒜, par exemple. Ce terme est moins daté que « discrimination » par exemple, et plus général que la constitution d'un « score » (qui est souvent une étape intermédiaire). Il ne doit pas être confondu avec la classification non-supervisée (comme la « classification ascendante hiérarchique ») qui est la constitution de classe homogène à partir d'une mesure de similarité (on utilisera parfois, dans ce cas, le terme de « constitution de classes », ou de « clusters »).], et dans le cas contraire, d'un problème de régression. Proches à l'arrivée, car dans les deux cas, on cherche à construire un « modèle », c'est à dire une fonction m:𝒳↦𝒴 qui sera interprétée comme une prévision.

Mais entre le départ et l'arrivée, il existe de réelles différences. Historiquement, les modèles économétriques s'appuient sur une théorie économique, avec le plus souvent des modèles paramétriques. On a alors recours aux outils classiques de l'inférence statistique (comme le maximum de vraisemblance, ou la méthode des moments) pour estimer les valeurs d'un vecteur de paramètres , dans un modèle paramétrique m_(·), par une valeur . Comme en statistique, avoir des estimateurs sans biais est important car on peut quantifier une borne inférieure pour la variance (borne de Cramér-Rao). La théorie asymptotique joue alors un rôle important (développements de Taylor, loi des grands nombres, et théorème central limite). En apprentissage statistique, en revanche, on construit souvent des modèles non-paramétriques, reposant presque exclusivement sur les données (i.e. sans hypothèse de distribution), et les méta-paramètres utilisés (profondeur de l'arbre, paramètre de pénalisation, etc) sont optimisés par validation croisée.

Au delà des fondements, si l'économétrie étudie abondamment les propriétés (souvent asymptotiques) de  (vu comme une variable aléatoire, grâce à la représentation stochastique sous-jacente), l'apprentissage statistique s'intéresse davantage aux propriétés du modèle optimal  m^⋆(·) (suivant un critère qui reste à définir), voire simplement m^⋆(x_i) pour quelques observations i jugées d'intérêt (par exemple dans une population de test). Le problème de choix de modèle est aussi vu sous un angle assez différent. Suivant la loi de Goodhart (« si une mesure devient un objectif, elle cesse d'être une mesure »), les économètres utilisent des critères de type AIC ou BIC pour choisir un modèle optimal (pénalisant la qualité d'ajustement d'un modèle par sa complexité, ex-post, lors de la phase de validation ou de choix), alors qu'en apprentissage statistique, c'est la fonction objectif qui tiendra compte d'une pénalisation, comme pour le lasso, ressemblant à une forme de pénalisation ex-ante.



 §.§ Applications


Avant de revenir sommairement sur l'évolution des modèles économétriques, c'est à Francis Galton que l'on doit le terme « régression », comme le rappelle KoenkerGalton. Si le terme est parfois devenu synonyme de « modèle é­co­no­mé­tri­que », il avait été introduit dans le contexte de « regression towards mediocraty inhereditary stature », pour reprendre le titre de l'article paru en 1886. Galton utilisait un modèle linéaire pour modéliser la taille moyenne d'un garçon (à l'âge adulte) en fonction de la taille de son père. Si cette technique de régression était connue par les économistes, il a fallu attendre les années 1930 pour voir surgir le concept de « modèle » économique. Comme le note Debreu, la première étape a été de formuler des affirmations économiques dans un language mathématique. Les différentes grandeurs sont vues comme des variables, et dans les années 1930, on verra apparaître les « statistical demand curves », pour reprendre la terminologie d'Henry Schultz. Cette approche statistique permettra d'aller plus loin que les travaux pionners de Engel qui étudiait empiriquement la relation entre la consommation et le revenu des ménages, par exemple, dans une approche uniquement descriptive.

Les modèles économétriques se sont développés en parallèle des modèles macro-économiques. Les premiers travaux de la Commission Cowles ont porté sur l'identification  des modèles économiques, et l'estimation de modèles à équations simultanées. Ces développements vont aboutir à un âge d'or de l'économétrie, dans les années 1960, où les modèles économétriques seront utilisés afin d'améliorer les prévisions macroéconomiques. On va alors voir apparaître tout un ensemble de « lois » qui sont souvent traduites comme des relations linéaires entre des grandeurs agrégées, telle que la « loi de Okun » introduite dans Okun qui postule une relation linéaire entre le variation du nombre de demandeurs d'emploi et de la croissance du PIB,

    ΔChômage_t=β_0+β_1Croissance_t+ε_t,

quand on étudie ces grandeurs au cours du temps (t), ou la loi de « Feldstein-Horioka » introduite dans FH qui suppose une relation linéaire entre les taux d'investissement et d'épargne, relativement au revenu national,

    investissement_i/revenu national_i=β_0+β_1
    épargne_i/revenu national_i+ε_i

quand on modèlise les liens entre les allocations investissement-épargne pour plusieurs pays (i). Cet âge d'or correspond aussi à un questionnement profond, suite à la critique de Lucas, s'interrogeant sur l'inefficacité de ces outils à expliquer et à prévoir des crises. La principale explication était alors le manque de fondement micro-économiques de ces modèles, ce qui donnera un second souffle aux modèles micro-économétriques. On pourra rappeler que cette critique dite « de Lucas » avait été formulée dans Orcutt, qui avançait l'idée que les données macroéconomiques posaient des problèmes insolubles d'identification. La solution passait forcément par de l'économétrie sur données individuelles (au lieu de données obtenues par aggrégation), ce qui sera reformulé quelques années plus tard par Koopmans. 

Malheureusement, les modèles micro-économétriques sont généralement plus complexes, car ils se doivent de tenir compte d'une éventuelle censure dans les données, avec par exemple le modèle introduit par Tobin, d'erreurs sur les variables (qui pourront être corrigées par des instruments avec une technique initiée par Reiersol) ou avoir été collectée avec un biais de sélection, avec les techniques proposées par Heckman. On notera que les économètres se sont beaucoup interrogés sur la manière dont les données étaient construites, et ne se sont jamais contentés de « construire des modèles ». Un exemple peut être l'évaluation des politiques publiques, largement détaillé dans Givord. Dans ce cas, en effet, deux écoles se sont opposées (initiant un débat que l'on verra resurgir tout au long de l'article sur les méthodes d'apprentissage statistique). La première, dite « structuraliste », cherchera à construire un modèle complet afin de décrire le comportement des agents é­co­no­mi­ques. La seconde, souvent qualifiée d'« empiriste », vise à tester l'effet d'une mesure sans pour autant expliciter les mécanismes sous-jacents. C'est ce qu'explique AngristKrueger, en affirmant « research
in a structuralist style relies heavily on economic theory to guide
empirical work [⋯] An alternative to structural modeling, [⋯] the
`experimentalist' approach, [⋯] puts front and center the problem
of identifying causal effects from specific events or situations ».

On peut aussi souligner que si l'approche de la Commission Cowles était très exigeante, en supposant le modèle connu, toute une littérature s'est développée en allégeant cette hypothèse, soit en essayant de choisir le bon modèle (avec les travaux de HendryK par exemple) ou en proposant de faire des moyennes de modèles (comme développé récemment par Lili). Et plus généralement, alors que l'analyse économétrique (en particulier à des fins de politique économique) s'est développée plus récemment autour de l'inférence causale, les techniques d'apprentissage machine ont été vues, traditionnellement, autour de la prédiction (où la recherche de corrélations suffisamment fortes entre variables suffit) d'où leur popularité dans des usages plus industriels de classification, comme la reconnaissance de caractères, de signature, d'images, ou de traduction, comme le rappelle Bishop. En biologie, ces techniques ont été appliquées pour créer des classifications d'espèces animales en fonction d'analyse d'ADN, ou dans le domaine militaire et sécuritaire pour l'identification de cibles ou de terroristes (potentiels). Il faudra attendre les années 1990 pour voir des applications en finance avec Altman par exemple, ou Herbrich pour une revue de littérature sur les applications potentielles en économie. Si des applications sont aujourd'hui nombreuses, et si ces techniques concurrencent les modèles de micro-économétrie (on pourra penser au scoring bancaire, à la détection de fraude fiscale ou assurantielle, ou à l'identification de prospects en marketing), les algorithmes d'apprentissage sont devenus très populaires en reconnaissance de parole, puis d'images, et plus récemment avec les applications en ligne et les applications aux jeux (d'échec, et plus récemment de go). Si l'économétrie s'est développée au confluent des mathématiques et de l'économie, l'apprentissage machine (que l'on pourrait avoir tendance à rapprocher de l'intelligence artificelle) s'est développé à la frontière des mathématiques et de l'informatique (avec des résultats fondamentaux en optimisation - en particulier autour des méthodes de gradient stochastique - et sur les espaces « sparse » ou « parcimonieux »).



 §.§ De la grande dimension aux données massives


Dans cet article, une variable sera un vecteur de ℝ^n, de telle sorte qu'en concaténant les variables ensemble, on puisse les stocker dans une matrice , de taille n× p, avec n et p potentiellement grands[Là encore, des extensions sont possibles, en particulier dans les données médicales avec des images de type IRM comme variables prédictives, ou des données climatiques avec des cartes en variables prédictives, ou plus généralement une variable tensorielle en dimension plus ou moins grande. Comme le montre Kolda il est toutefois possible de se ramener dans le cas usuel (de données sous formes de vecteurs) en utilisant la décomposition de Tucker.]. Le fait que n soit grand n'est, a priori, pas un problème en soi, au contraire. De nombreux théorèmes en économétrie et en statistique sont obtenus lorsque n→∞ (c'est la théorie asymptotique). En revanche, le fait que p soit grand est problématique, en particulier si p>n. Les deux dimensions sont à distinguer, car elles vont engendrer des problèmes relativement différents.

Portnoy a montré que l'estimateur du maximum de vraisemblance conserve la propriété de normalité asymptotique si p reste petit devant n, ou plus précisément, si p^2/n→ 0 lorsque n,p→∞. Aussi, il n'est pas rare de parler de grande dimension dès lors que p>√(n). Un autre concept important est celui de sparcité, qui repose non pas sur la dimension p mais sur la dimension effective, autrement dit le nombre de variables effectivement significatives. Il est alors possible d'avoir p>n tout en ayant des estimateurs convergents.

La grande dimension en terme de nombre de variables, p, peut faire peur à cause de la malédiction de la dimension, introduit par Bellman. L'ex­pli­ca­tion de cette malédiction est que le volume de la sphère unité, en dimension p, tend vers 0 lorsque p→∞. On dit alors que l'espace est « parcimonieux » - c'est à dire que la probabilité de trouver un point proche d'un autre devient de plus en plus faible (on pourrait parler d'espace « clairsemé »). Ou de manière duale, pour reprendre la formulation de HastieEtal, le volume qu'il convient de considérer pour avoir une proportion donnée d'observations augmente avec p. L'idée de réduire la dimension en considérant une analyse en composante principale peut paraître séduisante, mais l'analyse souffre d'un certain nombre de défauts en grande dimension. La solution est alors souvent la sélection de variables, qui pose le problème des tests multiples, ou du temps de calcul, pour sélectionner k variables parmi p, lorsque p est grand.

Pour reprendre la terminologie de Buhlmann, les problèmes que nous évoquons ici correspondent à ceux observés en grande dimension, qui est un problème essentiellement statistique. D'un point de vue informatique, on peut aller un peu plus loin, avec des données réellement massives (qui occupent énormément de place en mémoire). Dans ce qui précède, les données étaient stockées dans une matrice X, de taille n× p. Si cet objet formel est toujours bien défini, il peut y avoir des soucis à stocker une telle matrice, voire manipuler une matrice abondament utilisée en économétrie, X^ TX (matrice n× n). La condition du premier ordre (dans le modèle linéaire) est associée à la résolution de X^ T(Xβ-y)=0. En dimension raisonnable, on utilise la décomposition QR (c'est à dire la décomposition de Gram-Schmidt). En grande dimension, on peut utiliser des méthodes numériques de descente de gradient, où le gradient est approché sur un sous-échantillon de données (comme décrit par exemple dans Zinkevichetal) Cet aspect informatique est souvent oublié alors qu'il a été à la base de bon nombre d'avancées méthodologiques, en économétrie. Par exemple, HoerlKennard reviennent sur l'origine de l'utilisation de la régression Ridge: « Nous facturions  90$ par jour pour notre temps, mais avons dû facturer 450$ par heure d'ordinateur sur un Univac I (⋯) Avec cette machine, il a fallu 75 minutes de traitement pour inverser une matrice  40 × 40 en passant par une partition 4 × 4 de sous-matrices 10× 10, en utilisant des bandes magnétiques pour le stockage temporaire. Nous avons noté que les coefficients d'un régression linéaire calculés en utilisant les moindres carrés n'avaient pas toujours de sens. Les coefficients avaient tendance à être trop grands en valeur absolue, certains avaient même le mauvais signe, et ils pouvaient être instables avec de très petits changements dans les données (⋯) Comme la méthode que nous proposions attaquait l'une des vaches sacrées de la régression linéaire - les moindres carrés - nous avons fait face à une une résistance considérable». 






 §.§ Statistique computationnelle et non-paramétrique


L'objet de ce papier est d'expliquer les différences majeures entre l'économétrie et l'apprentissage statistique, correspondant aux deux cultures mentionnées par Breiman, lorsqu'il évoque en modélisation statistique la « data modeling culture » (reposant sur un modèle stochastique, comme la régression logistique ou le modèle de Cox) et la « algorithmic modeling culture » (reposant sur la mise en œuvre d'un algorithme, comme dans les forêts aléatoires ou les supports vecteurs machines, une liste exhaustive est présenté dans Shalev-Shwartz). Mais la frontière entre les deux est très poreuse. À l'intersection se retrouve, par exemple, l'économétrie non-paramétrique. Cette dernière repose sur un modèle probabiliste (comme l'é­co­no­mé­trie), tout en insistant davantage sur les algorithmes (et leurs performances) plutôt que sur des théorèmes asymptotiques.

L'économétrie non-paramétrique repose sur des décompositions dans des bases fonctionnelles. L'éco­no­mé­trie linéaire consiste à approcher la fonction m:↦[Y|=] par une fonction linéaire. Mais plus généralement, on peut considérer une décomposition dans une base fonctionnelle, et s'intéresser à une approximation obtenue sur un nombre fini de termes :

    m()=∑_j=0^∞ω_j g_j()     et   m()=∑_j=0^h^⋆ω_j g_j(),

où les poids ω_j sont estimés, alors que le nombre de composantes h^⋆ est optimisé. On retrouvera ici les modèles additifs (dits gam), par exemple, étudiés dans  Trevor. Une autre solution consiste à considérer un modèle simple, mais local. Par exemple un modèle constant, au voisinage de , obtenu en considérant seulement les observations proches de  :

    g()=∑_i=1^n ω_ y_i   par exemple   g()=1/n_∑_i:‖_i-‖≤ h y_i

où n_ est le nombre d'observations au voisinage de . En mettant des poids fonctions de la distance à , on retrouve ici le modèle obtenu par Nadaraya et Watson, ou les méthodes de régression locale.

Les différentes méthodes reposent sur des méta-paramètres - correspondant paramètres de lissage - c'est à dire h dans les exemples précédents. Pour un économètre, le paramètre « optimal » pour h est obtenu soit à l'aide de théorèmes asymptotiques, soit à l'aide de techniques de validation, comme en machine learning.  On obtient alors une valeur numérique, mais on n'a pas d'interprétation en lien avec la taille de l'échantillon, ou les variances des différentes grandeurs. Si les économistes ont toujours la culture du tableau présentant la« sortie de régression », les méthodes non-paramétriques sont utiles pour détecter des mauvaises spécifications, des non-prises en compte de nonlinéarité, ou d'effets croisées (et les outils de « machine learning » que nous allons voir peuvent probablement jouer le même rôle).



 §.§ Plan de l'article


Pour reprendre le titre de Varian, l'objet de cet article est de présenter les différences fondamentales entre l'économétrie et l'apprentissage machine, et surtout de voir comment ces deux techniques peuvent apprendre l'une de l'autre, dans un contexte où les bases de données deviennent massives. La Section <ref> reviendra sur la construction du modèle linéaire. Le modèle sera introduit ici à partir du modèle Gaussien « homoscédastique ». Ce modèle présente l'avantage d'avoir une élégante interprétation géométrique, en terme de projection sur le sous-espace des combinaisons linéaires des variables explicatives. La première extension que nous verrons est le passage du modèle linéaire à un modèle non-linéaire, tout en construisant un prédicteur linéaire. La seconde extension proposera de construire un modèle non-gaussien, pour modéliser une variable indicatrice ou un comptage Poissonnien, par exemple, donnant naissance aux modèles linéaires généralisés (construits pour des variables dans la famille exponentielle). 

Une fois rappelé l'origine des outils économétriques standards, dans la Section <ref> nous présenterons les outils et techniques développés dans le contexte de l'apprentissage machine. Si l'outil central des modèles économétriques est la distribution de la variable dépendante, Y, les techniques d'apprentissage reposent sur une fonction de perte, ℓ, représentant une « distance » entre la variable d'intérêt y, et le modèle m(·). Nous présenterons tout d'abord l'algorithme de boosting, reposant sur l'idée d'un apprentissage lent, en modélisant séquentiellement les résidus. Le danger des méthodes d'ap­pren­tis­sage est qu'il est aisé de construire un modèle « parfait », dont le pouvoir prédictif serait faible. Nous évoquerons alors les techniques de pénalisation, utilisées pour éviter le sur-apprentissage. Nous évoquerons en particulier les notions d'in-sample et out-of-sample, et les techniques de validation croisée. Pour conclure cette section, nous reviendrons sur les interprétations probabilistes des outils d'apprentissage, qui permettront de faire le lien entre les différentes approches, tout en restant sur une discussion générale sur la philosophie de ces deux cultures.

Après cette section sur la philosophie des méthodes de machine learning, nous reviendrons dans la section <ref> sur quelques algorithmes importants : les réseaux de neurones, les supports vecteurs machine (SVM) et enfin les méthodes de type arbres et forêts.

La Section <ref> proposera des exemples concrets de comparaison entre les différentes techniques, dans le cas de classifications (binaires) pour des variables y∈{0,1} (achat d'assurance, non-remboursement d'un crédit) et dans un contexte de régression (lorsque la variable d'intérêt n'est plus qualitative - ce que nous simplifierons en notant y∈ℝ). Nous reviendrons avant sur les courbes ROC, outils importants pour juger de la qualité d'un classifieur, malheureusement peu utilisés en économétrie.


Nous verrons en particulier les méthodes de bagging, forêts aléatoires ou boosting. Nous reviendrons aussi sur les méthodes de choix de modèles et des méta-paramètres. À travers ces exemples d'application, nous verrons comment les modèles de type machine learning peuvent être utilisés pour mieux détecter la mauvaise spécification des modèles de régression paramétriques, à cause de non-linéarités, et/ou d'intéractions manquées. 






§ ÉCONOMÉTRIE ET MODÈLE PROBABILISTE




L'importance des modèles probabilistes en économie trouve sa source dans les questionnements de Working et les tentatives de réponses apportées dans les deux tomes de Tinbergen. Ces derniers ont engendré par la suite énormément de travaux, comme le rappelle Duo dans son ouvrage sur les fondements de l'économétrie, et plus particulièrement dans le premier chapitre « The Probability Foundations of Econometrics ». Rappelons que Trygve Haavelmo a reçu le prix Nobel d’économie en 1989 pour sa « clarification des fondations de la théorie probabiliste de l'économétrie ». Car comme l'a montré Haavelmo (initiant un changement profond dans la théorie économétrique dans les années 1930, comme le rappelle le chapitre 8 de Morgan) l'économétrie repose fondamentalement sur un modèle probabiliste, et ceci pour deux raisons essentielles. Premièrement, l'utilisation de grandeurs (ou « mesures » ) statistiques telles que les moyennes, les erreurs-types et les coefficients de cor­ré­la­tion à des fins inférentielles ne peut se justifier que si le processus générant les données peut être exprimé en termes de modèle probabiliste. Deuxièmement, l'approche par les probabilitées est relativement générale, et se trouve être particulièrement adaptée à l'analyse des observations « dépendantes » et « non homogènes », telles qu'on les trouve souvent sur des données économiques. On va alors supposer qu'il existe un espace probabiliste (Ω,ℱ,ℙ) tel que les observations (y_i,_i) sont vues comme des réalisations de variables aléatoires (Y_i,_i). En pratique, la loi jointe du couple (Y,) nous intéresse toutefois peu : la loi de  est inconnue, et c'est la loi de Y conditionnelle à  qui nous intéressera. Dans la suite, nous noterons x une observation,  un vecteur d'observations,  X une variable aléatoire, et  un vecteur aléatoire et, abusivement,  pourra aussi désigner la matrice des observations individuelles (les _i), suivant le contexte.



 §.§ Fondements de la statistique mathématique


Comme le rappelle l'introduction de Vapnik98, l'inférence en statistique paramétrique est basée sur la croyance suivante:
le statisticien connaît bien le problème à analyser, en particulier, il connaît la loi physique qui génère les propriétés stochastiques des données, et la fonction à trouver s'écrit via un nombre fini de paramètres[On peut rapprocher cette approche de l'économétrie structurelle, telle que présentée par exemple dans Keen.]. Pour trouver ces paramètres, on adopte la méthode du maximum de vraisemblance. Le but de la théorie est de justifier cette approche (en découvrant et en décrivant ses propriétés favorables). On verra qu'en apprentissage, la philosophie est très différente, puisqu'on ne dispose pas d'informations a priori fiables sur la loi statistique sous-jacente au problème, ni-même sur la fonction que l'on voudrait approcher (on va alors proposer des méthodes pour construire une approximation à partir de données à notre disposition, pour reprendre Vapnik98).
Un « âge d'or » de l'inférence paramétrique, de 1930 à 1960, a posé les bases de la statistique mathématique, que l'on retrouve dans tous les manuels de statistique, y compris aujourd'hui. Comme le dit Vapnik98, le paradigme paramétrique classique est basé sur les trois croyances suivantes:



  * Pour trouver une relation fonctionnelle à partir des données, le statisticien est capable de définir un ensemble de fonctions, linéaires dans leurs paramètres, qui contiennent une bonne approximation de la fonction souhaitée. Le nombre de paramètres décrivant cet ensemble est petit.

  * La loi statistique sous-jacente à la composante stochastique de la plupart des problèmes de la vie réelle est la loi normale.
Cette croyance a été soutenue en se référant au théorème de limite centrale, qui stipule que dans de larges conditions la somme d'un grand nombre de variables aléatoires est approximée par la loi normale. 

  * La méthode du maximum de vraisemblance est un bon outil pour estimer les paramètres.

Nous reviendrons dans cette partie sur la construction du paradigme économétrique, directement inspiré de celui de la statistique inférentielle classique.




 §.§ Lois conditionnelles et vraisemblance


L'économétrie linéaire a été construite sous l'hypothèse de données individuelles, ce qui revient à supposer les variables  (Y_i,_i) in­dé­pen­dan­tes (s'il est possible d'imaginer des observations temporelles - on aurait alors un processus (Y_t,_t) - mais nous n'aborderons pas les séries temporelles dans cet article). Plus précisément, on va supposer que conditionnellement aux variables explicatives _i, les variables Y_i sont indépendantes. On va également supposer que ces lois conditionnelles restent dans la même famille paramétrique, mais que le paramètre est une fonction de . Dans le modèle linéaire Gaussien on suppose que :

    (Y|=)ℒ∼𝒩(μ(),σ^2)   avec   μ()=β_0+, et ∈ℝ^p.

On parle de modèle linéaire car [Y|=]=β_0+ est une combinaison linéaire des variables explicatives. C'est un modèle homoscédastique si [Y|=]=σ^2, où σ^2 est une constante positive. Pour estimer les paramètres, l'approche classique consiste à utiliser l'estimateur du Maximum de Vraisemblance, comme l'avait suggéré initialement Ronald Fisher. Dans le cas du modèle linéaire Gaussien, la log-vraisemblance s'écrit :

    logℒ(β_0,,σ^2|,) = -n/2log[2πσ^2] - 1/2σ^2∑_i=1^n (y_i-β_0-_i)^2.

Notons que le terme de droite, mesurant une distance entre les données et le modèle, va s'interpréter comme la déviance, dans les modèles linéaires généralisés.
On va alors poser :

    (β_0,,σ^2)=argmax{logℒ(β_0,,σ^2|,)}.

L'estimateur du maximum de vraisemblance est obtenu par minimisation de la somme des carrés des erreurs (estimateur dit des « moindres carrés » ) que nous retrouverons dans l'approche par machine learning.

Les conditions du premier ordre permettent de retrouver les équations normales, dont l'écriture matricielle est [-]=0, que l'on peut aussi écrire ()=. Si la matrice  est de plein rang colonne, alors on retrouve l'estimateur classique :

    =()= + ()

en utilisant une écriture basée sur les résidus (comme souvent en économétrie), y=+ε.
Le théorème de Gauss Markov assure que cette estimateur est l'estimateur linéaire sans biais de variance minimale. On peut alors montrer que ℒ∼𝒩(,σ^2 []), et en particulier :

    []=   et   []=σ^2 [].


En fait, l'hypothèse de normalité permet de faire un lien avec la statistique mathématique, mais il est possible de construire cet estimateur donné par l'équation (<ref>). Si on suppose que Y|=ℒ∼+ε, avec 𝔼[ε]=0, Var[ε]=σ^2 , Cov[ε,X_j]=0 pour tout j, alors  est un estimateur sans biais de  ([]=) et de variance minimale parmi les estimateurs sans biais linéaires, avec []=σ^2 []. De plus, cet estimateur est asymptotiquement normal  

    √(n)(-)ℒ→𝒩(0,Σ)    lorsque    n→∞



La condition d'avoir  une matrice  de plein rang peut être (numériquement) forte en grande dimension. Si elle n'est pas vérifiée, 
=() n'existe pas. Si 𝕀 désigne la matrice identité, notons toutefois que 
(+λ𝕀) existe toujours, pour λ>0. Cet estimateur est appelé l'estimateur Ridge de niveau λ (introduit dans les années 60 par Hoerl, et associé à une régularisation étudiée par Tikhonov). Cette estimateur apparaît naturellement dans un contexte d'économétrie Bayesienne (nous le reverrons dans la section suivante présentant les techniques de machine learning). 



 §.§ Les résidus


Il n'est pas rare d'introduire le modèle linéaire à partir de la loi des résidus, comme nous l'avions mentionné auparavant. Aussi, l'équation (<ref>) s'écrit aussi souvent :

    y_i=β_0+_i+ε_i

où les ε_i sont des réalisations de variables aléatoires i.i.d., de loi 𝒩(0,σ^2). On notera parfois ℒ∼𝒩(0,σ^2𝕀), sous une forme vectorielle. Les résidus estimés sont définis par :

    ε_i=y_i - [β_0+_i]

Ces résidus sont l'outil de base pour diagnostiquer la pertinence du modèle.

Une extension du modèle décrit par l'équation (<ref>) a été proposé pour tenir compte d'un éventuel caractère hétéroscédastique :

    (Y|=)ℒ∼𝒩(μ(),σ^2())

où σ^2() est une fonction positive des variables explicatives.
On peut réécrire ce modèle en posant :

    y_i=β_0+_i+σ^2(_i)·ε_i

où les résidus sont toujours i.i.d., mais de variance unitaire, 

    ε_i=y_i-[β_0+_i]/σ(_i).

Si l'écriture à l'aide des résidus est populaire en économétrie linéaire (lorsque la variable dépendante est continue), elle ne l'est toutefois plus dans les modèles de comptage, ou la régression logistique.

L'écriture à l'aide d'un terme d'erreur (comme dans l'équation (<ref>)) pose toutefois de nombreuses questions quant à la représentation d'une relation économique entre deux grandeurs. Par exemple, on peut supposer qu'il existe une relation (linéaire pour commencer) entre les quantités d'un bien échangé, q et son prix p. On peut ainsi imaginer une équation d'offre

    q_i=β_0+β_1 p_i+u_i

(u_i désignant un terme d'erreur) où la quantité vendue dépend du prix, mais de manière tout aussi légitime, on peut imaginer que le prix dépend de la quantité produite (ce qu'on pourrait appeler une équation de demande),

    p_i=α_0+α_1 q_i+v_i

(v_i désignant un autre terme d'erreur). 
Historiquement, le terme d'erreur dans l'équation (<ref>) a pu être inteprété comme une erreur idiosyncratique sur la variable y, les variables dites explicatives étant supposées fixées, mais cette interprétation rend souvent le lien entre une relation économique et un modèle économique compliqué, la théorie économique parlant de manière abstraite d'une relation entre grandeur, la modélisation économétrique imposant une forme spécifique (quelle grandeur est y et quelle grandeur est x) comme le montre plus en détails le chapitre 7 de Morgan.
 




 §.§ Géométrie du modèle linéaire gaussien




Définissons le produit scalaire dans ℝ^n, ⟨,⟩=, et notons ‖·‖ la norme euclidienne associée, ‖‖ = √() (notée ·_ℓ_2 dans la suite). Notons ℰ_ l'espace engendré par l'ensemble des combinaisons linéaires des composantes  (en rajoutant la constante). Si les variables explicatives sont linéairement indépendantes,  est de plein rang colonne et ℰ_ est un sous-espace de dimension p+1 de ℝ^n. Supposons à partir de maintenant que les variables  et la variable y sont ici centrées.
Notons qu'aucune hypothèse de loi n'est faite dans cette section, les propriétés géométriques découlent des propriétés de l'espérance et de la variance dans l'espace des variables de variance finie.

Avec cette notation, notons que le modèle linéaire s'écrit m()=⟨,β⟩. L'espace
ℋ_z={∈ℝ^k:m()=z} est un hyperplan (affine) qui sépare l'espace en deux.
Définissons l'opérateur de projection orthogonale sur ℰ_𝒳=ℋ_0, Π_𝒳=[]. Aussi, la prévision que l'on peut faire pour  est :

    =[]_
    Π_𝒳=Π_𝒳.

Comme =-=(𝕀-Π_𝒳)=Π_𝒳^⊥, on note que ⊥, que l'on interprétera en disant que les résidus sont un terme d'innovation, imprévisible, au sens où Π_𝒳=0.

Le théorème de Pythagore s'écrit ici :

    ‖‖^2=‖Π_𝒳‖^2+‖Π_𝒳^⊥‖^2=‖Π_𝒳‖^2+‖-Π_𝒳‖^2=‖‖^2+‖ε‖^2

qui se traduit classiquement en terme de somme de carrés :

    ∑_i=1^n y_i^2_n×variance totale=∑_i=1^n y_i^2_n×variance expliquée
    +
    ∑_i=1^n (y_i-y_i)^2_n×variance résiduelle

Le coefficient de détermination, R^2 (nous reviendrons sur ce coefficient dans la section <ref>) s'interprête alors comme le carré du cosinus de l'angle θ entre  et Π_𝒳 :

    R^2=‖Π_𝒳‖^2/‖‖^2=1-‖Π_𝒳^⊥‖^2/‖‖^2=cos^2(θ).



Une application importante a été obtenue par FrishWaugh, lorsque l'on partitionne les variables explicatives en deux groupes, =[ _1 |_2], de telle sorte que la régression devient :

    =β_0+ _1_1+_2_2+ε

FrishWaugh ont montré qu'on pouvait considérer deux projections successives. En effet, si _2^⋆=Π_𝒳_1^⊥ et _2^⋆=Π_𝒳_1^⊥_2, on peut montrer que

    _2=[_2^⋆_2^⋆]_2^⋆_2^⋆

Autrement dit, l'estimation globale est équivalente à l'estimation indépendante des deux modèles si _2^⋆=_2, c'est à dire _2∈ℰ__1^⊥, que l'on peut noter _1⊥_2. On obtient ici le théorème de Frisch-Waugh qui garantie que si les variables explicatives entre les deux groupes sont orthogonales, alors l'estimation globale est équivalente à deux régressions indépendantes, sur chacun des jeux de variables explicatives. Ce qui est un théorème de double projection, sur des espaces orthogonaux. Beaucoup de résultats et d'interprétations sont obtenus par des interprétations géométriques (liées fondamentalement aux liens entre l'espérance conditionnelle et la projection orthogonale dans l'espace des variables de variance finie).

Cette vision géométrique permet de mieux comprendre le problème de la sous-identification, c'est à dire le cas où le vrai modèle serait y_i=β_0+_1_1+_2_2+ε_i, mais le modèle estimé est y_i=β_0+_1b_1+η_i. L'estimateur du maximum de vraisemblance de b_1 est :

    b_1
       =   (X_1^TX_1)^-1X_1^T
       =   (X_1^TX_1)^-1X_1^T
    [X_1,iβ_1
        + X_2,iβ_2 + ε] 
       =    (X_1^TX_1)^-1X_1^TX_1β_1
    + (X_1^TX_1)^-1X_1^TX_2β_2
    + (X_1^TX_1)^-1X_1^Tε
       =   β_1 +  (X_1'X_1)^-1X_1^TX_2β_2_β_12
    +(X_1^TX_1)^-1X_1^Tε_ν_i

de telle sorte que [b_1]=_1+_12, le biais étant nul uniquement dans le cas où X_1^TX_2=0 (c'est à dire X_1 ⊥X_2): on retrouve ici une conséquence du théorème de Frisch-Waugh.

En revanche, la sur-identification correspond au cas où le vrai modèle serait y_i=β_0+_1_1+ε_i, mais le modèle estimé est y_i=β_0+_1b_1+_2b_2+η_i.
Dans ce cas, l'estimation est sans biais, au sens où 𝔼(b_1)=β_1 mais l'estimateur n'est pas efficient. Et comme nous l'avons vu dans la section précédente, il n'est pas rare d'avoir des valeurs de b_2 qui sont considérées comme significativement non-nulles. Nous évoquerons dans la section suivante une méthode efficace de choix de variables (et éviter la sur-identification).



 §.§ Du paramétrique au  non-paramétrique


La réécriture de l'équation (<ref>) sous la forme

    = 
    =[]_Π_𝒳
 
permet de voir la prévision directement comme une transformation linéaire des observations. De manière plus générale, on peut obtenir un prédicteur linéaire en considérant m()=_, où _ est un vecteur de poids, qui dépendent de , interprété comme un vecteur de lissage. En utilisant les vecteurs __i, calculés à partir des _i, on obtient une matrice  S de taille n× n, et =S.
Dans le cas de la régression linéaire décrite auparavant, _=[], et classiquement, trace( S) est le nombre de colonnes de la matrice  (le nombre de variables explicatives). Dans ce contexte de prédicteurs linéaires,  trace( S) est souvent vu comme un équivalent au nombre de paramètres (ou complexité, ou dimension, du modèle), et ν=n-trace( S) est alors le nombre de degrés de liberté (comme défini dans Ruppert et Simonoff). Le principe de parcimonie[« pluralitas non est
ponenda sine necessitate » pour reprendre le principe énoncé par Guillaume d'Occam (les multiples ne
doivent pas être utilisés sans nécessité).] consiste à minimiser cette dimension (la trace de la matrice S) autant que faire se peut. Mais dans le cas général, cette dimension est plus complexe à définir. Notons que l'estimateur introduit par Nadaraya et Watson, dans le cas d'une régression non-paramétrique simple, s'écrit également sous cette forme puisque

    m_h(x)=_x=∑_i=1^n s_x,iy_i    avec    s_x,i=K_h(x-x_i)/K_h(x-x_1)+⋯+K_h(x-x_n),

où K(·) est une fonction noyau, qui attribue une valeur d'autant plus faible que x_i est proche de x, et h>0 est la fenêtre de lissage. 




L'introduction de ce méta-paramètre h pose un soucis, car il convient de le choisir judicieusement. En faisant des développement limités, on peut montrer que si X a pour densité f,

    biais[m_h(x)]=𝔼[m_h(x)]-m(x)∼h^2(
    C_1 /2m”(x)+C_2 m'(x)f'(x)/f(x))
     et Var[m_h(x)]∼C_3/nhσ(x)/f(x)

pour des constantes que l'on peut estimer (voir Simonoff par exemple). Ces deux fonctions évoluent inversement en fonction de h, comme le rappelle la Figure <ref>. L'idée naturelle est alors de chercher à minimiser l'erreur quadratique moyenne, le MSE, biais[m_h(x)]^2+Var[m_h(x)], ce qui donne une valeur optimale pour h de la forme h^⋆=O(n^-1/5), ce qui rappelle la règle de Silverman. En plus grande dimension, pour des variables  continues, on peut utiliser un noyau multivarié, de fenêtre matricielle H, 

    𝔼[m_H(x)]∼ m(x)
    +C_1/2trace(H^ Tm”(x)H)+C_2m'(x)^ THH^ T∇ f(x)/f(x) et Var[m_H(x)]∼C_3/n det(H)σ(x)/f(x).

Si H est une matrice diagonale, avec le même terme h sur la diagonale, alors h^⋆=O(n^-1/(4+dim(x)))]. Cela dit, en pratique, on sera davantage intéressé par la version intégrée de l'erreur quadratique,

    MISE(m_h)=𝔼[MSE(m_h(X))]=∫ MSE(m_h(x))dF(x),

dont on peut montrer que

    MISE[m_h]   ∼   h^4/4(∫ x^2k(x)dx)^2∫[m”(x)+2m'(x)f'(x)/f(x)]^2dx^biais^2 +σ^2/nh∫ k^2(x)dx ·∫dx/f(x)^variance,

lorsque n→∞ et nh→∞. On retrouve ici une relation asymptotique qui rappelle l'ordre de grandeur de Silverman,

    h^⋆ =n^-1/5(C_1∫dx/f(x)/C_2∫[m”(x)+2m'(x)f'(x)/f(x)]dx)^1/5,

sauf que beaucoup de termes ici sont inconnus. On verra, le machine-learning propose des techniques computationnelles, lorsque l'économètre avait pris l'habitude de chercher des propriétés asymptotiques.







 §.§ Famille exponentielle et modèles linéaires


Le modèle linéaire Gaussien est un cas particulier d'une vaste famille de modèles linéaires, obtenu lorsque la loi conditionnelle de Y appartient à la famille exponentielle

    f(y_i|θ_i,ϕ)=exp(y_iθ_i-b(θ_i)/a(ϕ)+c(y_i,ϕ))   avec   θ_i=ψ(x_i^Tβ).

Les fonctions a, b et c sont spécifiées en fonction du type de
loi exponentielle (étudiée abondamment en statistique depuis les Darmois, comme le rappelle Brown), et ψ est une fonction bijective que se donne l'utilisateur. 
La log-vraisemblance a alors une expression relative simple

    logℒ(θ,ϕ|y)=∏_i=1^n log f(y_i|θ_i,ϕ)
    =∑_i=1^ny_iθ_i-∑_i=1^nb(θ_i)/a(ϕ)+
    ∑_i=1^n c(y_i,ϕ)

et la condition du premier ordre s'écrit alors

    ∂logℒ(θ,ϕ|y)/∂β = X^TW^-1[y-]=0

pour reprendre les notations de Muller, où W est une matrice de poids (qui dépend de ). Compte tenu du lien entre θ et l'espérance de Y, au lieu de spécifier la fonction ψ(·), on aura plutôt tendance à specifier la fonction de lien g(·) définie par 

    y=m()=[Y|=]=g().

Pour la régression linéaire Gaussienne on prendra un lien Identité, alors que pour la régression de Poisson, le lien naturel (dit canonique) est le lien logarithmique. Ici, comme W dépend de  (avec  W=diag(∇ g()Var[])) il n'existe en général par de formule explicite pour l'estimateur du maximum de vraisemblance. Mais un algorithme itératif permet d'obtenir une approximation numérique. En posant

    =g()+(-)·∇ g()

correspondant au terme d'erreur d'un développement de Taylor à l'ordre 1 de g, on obtient un algorithme de la forme

    β_k+1 = [X^ TW_k^-1X]^-1X^ TW_k^-1z_k

En itérant, on notera β=β_∞, et on peut montrer que - moyennant quelques hypothèses techniques (cf Muller) - cet estimateur est asymptotiquement Gaussien, avec

    √(n)(β-β)ℒ→𝒩(0,I(β)^-1),

où numériquement I(β)=ϕ·  [X^ TW_∞^-1X].

D'un point de vue numérique toujours, on résout la condition du premier ordre, et la loi de Y n'intervient pas réellement. Par exemple, on peut estimer une « régression de Poisson » même lorsque y∈ℝ_+, pas nécessairement y∈ℕ. Autrement dit, la loi de Y n'est qu'une interprétation donnée ici, et l'algorithme pourrait être introduit de manière différente (comme nous le verrons dans la section suivante), sans forcément avoir de modèle probabiliste sous-jacent.




 §.§ Régression logistique


La régression logistique est le modèle linéaire généralisé obtenu avec une loi de Bernoulli, et une fonction de lien qui est la fonction quantile d'une loi logistique (ce qui correspond au lien canonique au sens de la famille exponentielle). Compte tenu de la forme de la loi de Bernoulli, l'économétrie propose un modèle pour y_i∈{0,1}, dans lequel le logarithme de la cote suit un modèle linéaire :

    log(ℙ[Y=1|=]/ℙ[Y≠ 1|=])=β_0+,

ou encore :

    [Y|=]=ℙ[Y=1|=]=e^β_0+/1+e^β_0+=H(β_0+),    où  H(·)=exp(·)/1+exp(·),

correspondant à la fonction de répartition de la loi logistique.
L'estimation de (β_0,) se fait par maximisation de la vraisemblance :

    ℒ=∏_i=1^n (e^x_i^ Tβ/1+e^x_i^ Tβ)^y_i(1/1+e^x_i^ Tβ)^1-y_i

On continuera à parler des modèles linéaires car les courbes d'isoprobabilités sont ici les hyperplans parallèles b_0+. À ce modèle, popularisé par Berkson1, certains préfèront le modèle probit (comme le raconte Berkson2), introduit par Bliss. Dans ce modèle :

    [Y|=]=ℙ[Y=1|=]=Φ(β_0+),

où Φ désigne la fonction de répartition de la loi normale centrée réduite. Ce modèle présente l'avantage d'avoir un lien direct avec le modèle linéaire Gaussien, puisque 

    y_i=1(y_i^⋆>0)   avec   
    y_i^⋆=β_0++ε_i

où les résidus sont Gaussiens, de loi 𝒩(0,σ^2). Un alternative est d'avoir des résidus centrés de variance unitaire, et de considérer une modélisation latente de la forme y_i=1(y_i^⋆>ξ) (où ξ sera à fixer). On le voit, ces techniques sont fondamentalement liées à un modèle stochastique sous-jacent. Mais dans la section <ref>, nous présenterons plusieurs techniques alternatives - tirées de la littérature en apprentissage - pour ce problème de classification (avec deux classes, ici 0 et 1).



 §.§ Régression en grande dimension


Comme nous l'avions mentionné auparavant, la condition du premier ordre X(Xβ-y)=0 se résout numériquement en effectuant une décomposition QR, pour un coût en O(np^2) opérations (où p est le rang de ). Numériquement, ce calcul peut être long (soit parce que p est grand, soit parce que n est grand), et une stratégie plus simple peut être de faire du sous-échantillonnage. Soit n_s≪ n, et considérons un sous-échantillon de taille n_s de {1,⋯,n}. Alors  β_s=(_s_s)^-1_sy_s est une bonne approximation de β comme le montre Dhilonetal2013. Cet algorithme est toutefois dangereux si certains points ont un pouvoir de levier important (i.e. L_i=_i()^-1_i). Tropp2011 propose de transformer les données (de manière linéaire), mais une approche plus populaire est de faire du sous-échantillonnage non uniforme, avec une probabilité liée à l'influence des observations (définie par I_i=ε_iL_i/(1-L_i)^2, et qui malheureusement ne peut être calculée qu'une fois le modèle estimé).


De manière générale, on parlera de données massives lorsque  la table de données de taille n × p ne tient pas en mémoire RAM de l'ordinateur. Cette situation est souvent rencontrée en apprentissage statistique de nos jours avec très souvent p ≪ n. C'est la raison pour laquelle, en pratique de nombreuses bibliothèques d'algorithmes assimilées à de l'apprentissage machine[comme, par exemple, celles du langage Python.] utilisent des méthodes itératives pour résoudre la condition du premier ordre. Lorsque le modèle paramétrique à calibrer est effectivement convexe et semi-différentiable, il est possible d'utiliser par exemple la méthode de descente de gradient stochastique comme le suggère bottou. Ce dernier permet de s'affranchir à chaque itération du calcul du gradient sur chaque observation de notre base d'apprentissage. 
Plutôt que d'effectuer une descente moyenne à chaque itération, on commence par tirer (sans remise) une observation X_i parmi les n disponibles. On corrige ensuite les paramètres du modèle de sorte à ce que la prédiction faite à partir de X_i soit la plus proche possible de la vraie valeur y_i. On réitère ensuite la méthode jusqu'à avoir parcourue l'ensemble des données. Dans cet algorithme il y a donc autant d'itération que d'observations. Contrairement à l'algorithme de descente de gradient (ou méthode de Newton) à chaque itération un seul vecteur de gradient est calculé (et non plus n). Il est néanmoins parfois nécessaire d'exécuter cette algorithme plusieurs fois pour augmenter la convergence des paramètres du modèle.

    
Si l'objectif est par exemple de minimiser l'erreur quadratique ℓ entre l'estimateur f_β(X) et y l'algorithme peut se résumer ainsi : 

Etape 0: Mélange des données

Etape d'itérations: Pour t=1,...,n, on tire i ∈{ 1,⋯, n } sans remise, on pose

    β^t+1 = β^t - γ_t∂l(y_i,f_β^t(X_i)) /∂β


  
Cet algorithme peut être réitérée plusieurs fois dans son ensemble selon le besoin de l'utilisateur.
L'avantage de cette méthode est qu'à chaque itération, il n'est pas nécessaire de calculer le gradient sur toutes les observations (plus de somme). Elle est donc adaptée aux bases de données volumineuses. Cet algorithme s'appuie sur une convergence en probabilité vers un voisinage de l'optimum (et non pas l'optimum lui même).













 §.§ Qualité d'un ajustement et choix de modèle


Dans le modèle linéaire Gaussien, le coefficient de détermination - noté R^2 - est souvent utilisé comme mesure de la qualité d'ajustement. Compte tenu de la formule de décomposition de la variance

    1/n∑_i=1^n (y_i-y̅)^2_variance totale=1/n∑_i=1^n (y_i-y_i)^2_variance résiduelle+1/n∑_i=1^n (y_i-y̅)^2_variance expliquée

on définit le R^2 comme le ratio de variance expliquée et de la variance totale, autre interpétation du coefficient que nous avions introduit à partir de la géométrie des moindres carrés.

    R^2 = ∑_i=1^n (y_i-y̅)^2-∑_i=1^n (y_i-y_i)^2/∑_i=1^n (y_i-y̅)^2

Les sommes des carrés d'erreurs dans cette écriture peut se réécrire comme une log-vraisemblance. Or rappelons qu'à une constante près, dans les modèles linéaires généralisés, la déviance est définie par

    Deviance() = -2log[ℒ]

que l'on peut aussi noter  Deviance(). On peut définir une déviance nulle comme celle obtenue sans utiliser les variables explicatives , de telle sorte que y_i=y. On peut alors définir, dans un contexte plus général

    R^2=
    Deviance(y)-Deviance(y)/Deviance(y)=1-
    Deviance(y)/Deviance(y).


Toutefois, cette mesure ne peut être utilisée pour choisir un modèle, si on souhaite avoir au final un modèle relativement simple, car elle augmente artificiellement avec l'ajout de variables explicatives sans effet significatif. On aura alors tendance à préférer le R^2 ajusté

    R̅^2 = 1-(1-R^2)n-1  n-p = R^2-(1-R^2)p-1  n-p_pénalisation,

où p est le nombre de paramètres du modèle (noté plus généralement ν dans la section <ref>). À la mesure de la qualité de l'ajustement, on va pénaliser les modèles trop complexes.

Cette idée va se retrouver dans le critère d'Akaike, où AIC = Deviance + 2· p ou dans le critère de Schwarz, BIC = Deviance + log(n)· p. En grande dimension (typiquement p>√(n)), on aura tendance à utiliser un AIC corrigé, défini par

    AICc = Deviance + 2· p·n/n-p-1


Ces critères sont utilisés dans les méthodes dites « stepwise », introduisant les méthodes ensemblistes. Dans la méthode dite « forward », on commence par régresser sur la constante, puis on ajoute une variable à la fois, en retenant celle qui fait le plus baisser le critère AIC, jusqu'à ce que rajouter une variable augmente le critère AIC du modèle. Dans la méthode dite « backward », on commence par régresser sur toutes les variables, puis on enlève une variable à la fois, en retirant celle qui fait le plus baisser le critère AIC, jusqu'à ce que retirer une variable augmente le critère AIC du modèle. 

















Une autre justification de cette notion de pénalisation (nous reviendrons sur cette idée en apprentissage) peut être la suivante. Considérons un estimateur dans la classe des prédicteurs linéaires,

    ℳ={ m: m()=s_h() oùS=(s(_1),⋯,s(_n) est la matrice de lissage}

et supposons que y=m_0(x)+ε, avec 𝔼[ε]=0 et Var[ε]=σ^2𝕀, de telle sorte que 
m_0(x)=𝔼[Y|X=x].
D'un point de vue théorique, le risque quadratique, associé à un modèle estimé m, s'écrit 

    ℛ(m)=𝔼[(Y-m(X))^2]
    =
    𝔼[(Y-m_0(X))^2]_erreur+
    𝔼[(m_0(X)-𝔼[m()])^2]_biais+
    𝔼[(𝔼[m()]-m())^2]_variance

si m_0 désigne le vrai modèle. Le premier terme est parfois appelé « erreur de Bayes », et ne dépend pas de l'estimateur retenu, m.

Le risque empirique quadratique, associé à un modèle m, est ici :

    ℛ_n(m)=1/n∑_i=1^n (y_i-m(x_i))^2 =1/ny-m(x)^2

(par convention). On reconnaît ici l'erreur quadratique moyenne, mse, qui donnera plus généralement le « risque » du modèle m quand on utilise une autre fonction de perte (comme nous le discuterons dans la partie suivante). Notons que:

    𝔼[ℛ_n(m)]=1/nm_0(x)-m(x)^2+1/n𝔼(y-m_0(x)^2)

On peut montrer que :

    n𝔼[ℛ_n(m)]
    =𝔼(
    y-m(x)^2)
    =
    (𝕀-S)m_0^2+σ^2
    𝕀-S^2

de telle sorte que le (vrai) risque de m est :

    ℛ_n(m)=
    𝔼[ℛ_n(m)]+2σ^2/ntrace(S).

Aussi, si trace(S)≥ 0, le risque empirique sous-estime le vrai risque de l'estimateur. On reconnaît ici le nombre de degrés de liberté du modèle, le terme de droite correspondant au C_p de Mallow, introduit dans Mallows utilisant non pas la déviance mais le R^2.



 §.§ Économétrie et tests statistiques


Le test le plus classique en économétrie est probablement le test de significativité, correspondant à la nullité d'un coefficient dans un modèle de régression linéaire. Formellement, il s'agit du test de H_0:β_k=0 contre H_1:β_k≠ 0.
Le test de Student, basé  sur la statistique t_k=β_k/se_β_k, permet a priori de trancher entre les deux alternatives, à l'aide de la p-value du test, définie par ℙ[|T|>|t_k|] avec Tℒ∼ Std_ν, où ν est le nombre de degrés de liberté du modèle (ν=p+1 pour le modèle linéaire standard).
En grande dimension, cette statistique est néanmoins d'un intérêt très limité, compte tenu d'un FDR (False Discovery Ratio) important. Classiquement, avec un niveau de significativité α=0.05, 5% des variables sont faussement significatives. Supposons que nous disposions de p=100 variables explicatives, mais que 5 (seulement) sont réellement significatives. On peut espérer que ces 5 variables passeront le test de Student, mais on peut aussi s'attendre à
ce que 5 variables supplémentaires (test faussement positif) ressortent. On aura alors 10 variables perçues comme significatives, alors que seulement la moitié le sont, soit un ratio FDR de 50%. Afin d'éviter cet écueil récurent dans les tests multiples, il est naturel d'utiliser la procédure de Benjamini.











 §.§ Quitter la corrélation pour quantifier un effet causal


Les modèles économétriques sont utilisés pour mettre en oeuvre des politiques publiques. Il est alors fondamental de bien comprendre les mécanismes sous-jacents pour savoir quelles variables permettent effectivement d'agir sur une variable d'intérêt. Mais on passe alors dans une autre dimension importante de l'économétrie. C'est à Jerry Neyman que l'on doit les premiers travaux sur l'identification de mécanismes causaux, c'est Rubin qui a formalisé le test, appelé « modèle causal de Rubin » dans Holland. Les premières approches autour de la notion de causalité en économétrie se sont faites avec l'utilisation des variables instrumentales, des modèles avec discontinuité de régression, l'analyse de différences dans les différences, ainsi que des expériences naturelles ou pas. La causalité est généralement déduite en comparant l'effet d'une politique - ou plus généralement d'un traitement - avec son contrefactuel, idéalement donné par un groupe témoin, aléatoire.  L'effet causal du traitement est alors défini comme Δ=y_1-y_0, c'est à dire la différence entre ce que serait la situation avec traitement (noté t=1) et sans traitement (noté t=0). Le souci est que seul y=t· y_1+(1-t)y_0 et t sont observés. Autrement dit l'effet causal de la variable t sur y n'est pas observé (puisque seule une des deux variables potentielles - y_0 ou y_1 est observée pour chaque individu), mais il est aussi individuel, et donc fonction de covariables . Généralement, en faisant des hypothèses sur la distribution du triplet (Y_0,Y_1,T), certains paramètres de la distribution de l'effet causal deviennent identifiables, à partir de la densité 
des variables observables (Y, T). Classiquement, on sera intéressé par les moments de cette distribution, en particulier l’effet moyen du traitement dans la population, 𝔼[Δ], voire juste l'effet moyen du traitement en cas de traitement 𝔼[Δ|T=1]. Si le résultat (Y_0,Y_1) est indépendant de la variable d’accès au traitement T, on peut montrer que
𝔼[Δ]=𝔼[Y|T=1]-𝔼[Y|T=0]. Mais si cette hypothèse d'indépendance n'est pas vérifiée, on a un biais de sélection, souvent associé à 𝔼[Y_0|T=1]-𝔼[Y_0|T=0]. RosenbaumRubin proposent d'utiliser un score de propension à être traité, p()=ℙ[T=1|X=x], en notant que si la variable Y_0 est indépendante de l'accès au traitement T conditionnellement aux variables explicatives , alors elle est indépendante de T conditionnellement au score p() : il suffit de les apparier à l'aide de leur score de propension.
Heckman2 propose ainsi un estimateur à noyau sur le score de propension, ce qui permet d'avoir simplement un estimateur de l'effet du traitement, conditionnellement au fait d’être traité.




§ PHILOSOPHIE DES MÉTHODES DE MACHINE LEARNING


Parallèlement à ces outils développés par et pour des économistes, toute une littérature a été développée sur des questions similaires, centrées autour de la prévision. Pour Breiman, une première différence vient du fait que la statistique s'est développée autour du principe d'inférence (ou d'expliciter la relation liant y aux variables ) alors qu'une autre culture s'intéresse avant tout à la prédiction. Dans une discussion qui suit l'article, David Cox l'affirme très clairement « predictive success (...) is not the primary basis for model choice ». 
Nous allons présenter les fondements des techniques du machine learning (les exemples d'algorithmes étant présentés dans les sections suivantes). Le point important, comme nous allons le voir, est que la principale préoccupation de l'apprentissage machine est liée aux propriétés de généralisation d'un modèle, c'est-à-dire sa performance - selon un critère choisi a priori - sur des données nouvelles, et donc des tests hors échantillon.



 §.§ Apprentissage par une machine


Aujourd'hui, on parle de « machine learning » pour décrire tout un ensemble de techniques, souvent computationnelles, alternatives à l'approche décrite auparavant (correspondant à l'économétrie classique). Avant de les caractériser autant que possible, notons juste qu'historiquement d'autres noms ont pu être donnés. Par exemple, Friedman97 propose de faire le lien entre la statistique (qui ressemble beaucoup aux techniques économétriques - test d'hypothèses, ANOVA, régression linéaire, logistique, GLM, etc) et ce qu'il appelait alors « data mining » (qui englobait alors les arbres de décisions, les méthodes des plus proches voisins, les réseaux de neurones, etc.). Le pont qu'il contribuera à construire correspond aux techniques d'apprentissages statistiques, décrites dans HastieEtal, mais l'apprentissage machine est un très vaste champ de recherche.

L'apprentissage dit « naturel » (par opposition à celui d'une machine) est celui des enfants, qui apprennent à parler, à lire, à jouer. Apprendre à parler signifie segmenter et catégoriser des sons, et les associer à des significations. Un enfant apprend aussi simultanément la structure de sa langue maternelle et acquiert un ensemble de mots décrivant le monde qui l'entoure. Plusieurs techniques sont possible, allant d'un apprentissage par coeur, par généralisation, par découverte, apprentissage plus ou moins supervisé ou autonome, etc. L'idée en intelligence artificielle est de s'inspirer du fonctionnement du cerveau pour apprendre, pour permettre un apprentissage « artificiel » ou « automatique », par une machine. Une première application a été d'apprendre à une machine à jouer à un jeux (tic-tac-toe, échecs, go, etc). Une étape indispensable est d'expliquer l'objectif qu'il doit atteindre pour gagner. Une approche historique a été de lui apprendre les règles du jeu. Si cela permet de jouer, cela ne permettra pas à la machine de bien  jouer. En supposant que la machine connaisse les règles du jeu, et qu'elle a le choix entre plusieurs dizaines de coups possible, lequel doit-elle choisir ? L'approche classique en intelligence artificielle utilise l'algorithme dit min-max utilisant une fonction d'évaluation : dans cet algorithme, la machine effectue une recherche en avant dans l'arbre des coups possibles, aussi loin que les ressources de
calcul le lui permettent (une dizaine de coups aux échecs, par exemple). Ensuite, elle calcule différents critères (qui lui ont été indiqué au préalable) pour toutes les positions (nombre de pièces prises, ou perdues, occupation du centre, etc. dans notre exemple du jeu d'échec), et finalement, la machine joue le coup qui lui permet de maximiser son gain. Un autre exemple peut être celui de la classification et de la reconnaissance d'images ou de formes. Par exemple, la machine doit identifier un chiffre dans une écriture manuscrite (chèque, code postal). Il s'agit de prédire la valeur d'une variable y, en sachant qu'a priori y∈{0,1,2,⋯,8,9}. Un stratégie classique est de fournir à la machine des bases d'apprentissage, autrement dit ici des millions d'images labélisées (identifiées) de chiffres manuscrits. Une stratégie simple et naturelle est d'utiliser un critère de décision basé sur les plus proches voisins dont on connaît l'étiquette (à l'aide d'une métrique prédéfinie).

La méthode des plus proches voisins (« k-nearest neighbors ) peut être décrit de la manière suivante : on considère (comme dans la partie précédante) un ensemble de n observations, c'est à dire des paires (y_i,x_i) avec x_i∈ℝ^p. Considérons une distance Δ sur ℝ^p (la distance Euclienne ou la distance de Mahalanobis, par exemple). Étant donnée une nouvelle observation x∈ℝ^p, supposons les observations ordonnées en fonction de la distance entre les x_i et x, au sens où

    Δ(x_1,x)≤Δ(x_2,x)≤⋯≤Δ(x_n,x)

alors on peut considérer comme prédiction pour y la moyenne des k plus proches voisins,

    m_k(x)=1/k∑_i=1^k y_i.

L'apprentissage fonctionne ici par induction, à partir d'un échantillon (appelé base d'apprentissage).

Le Machine Learning englobe ces algorithmes qui donnent aux ordinateurs la capacité d'apprendre sans être explicitement programmé (comme l'avait défini Arthur Samuel en 1959). La machine va alors explorer les données avec un objectif précis (comme chercher les plus proches voisins dans l'exemple que nous venons de décrire). Tom Mitchell a proposé une définition plus précise en 1998 : on dit qu'un programme d'ordinateur apprend de l'expérience E par rapport à une tâche T et une mesure de performance P, si sa performance sur T, mesurée par P, s'améliore avec l'expérience E. La tâche T peut être un score de défaut par exemple, et la performance P peut être le pourcentage d'erreurs commise. Le système apprend si le pourcentage de défauts prédit augmente avec l'expérience.

On le voit, l'apprentissage machine est fondamentalement un problème d'optimisation d'un critère à partir de données (dites d'apprentissage). Nombreux sont les ouvrages de programmation qui proposent des algorithmes, sans jamais faire mention d'un quelconque modèle probabiliste. Dans Watt par exemple, il n'est fait mention du mot « probabilité » qu'une seule fois, avec cette note de bas de page qui surprendra et fera sourire les économètres, « logistic regression can also be interpreted from a probabilistic perspective » (page 86). Mais beaucoup d'ouvrages récents proposent une relecture des approches d'apprentissage machine à l'aide de théories probabilistes, suite aux travaux de Vaillant et Vapnik. En proposant le paradigme de l'apprentissage « probablement à peu près correct » (PAC), une saveur probabiliste a été rajouté à l'approche jusqu'alors très computationnelle, en quantifiant l'erreur de l'algorithme d'apprentissage (dans un problème de classification).



 §.§ Le tournant des années 80/90 et le formalisme probabiliste


On dispose d'un échantillon d'apprentissage, avec des observations (x_i,y_i) où les variables y sont dans un ensemble 𝒴. Dans le cas de la classification, 𝒴={-1,+1}, mais on peut imaginer un ensemble relativement général. Un prédicteur est une fonction m à valeurs dans 𝒴, permettant d'étiqueter (ou de classer) les nouvelles observations à venir. On suppose que les étiquettes sont produites par un classifieur f appelé cible. Pour un statisticien, cette fonction serait le vrai modèle. Naturellement, on veut construire m le plus proche possible de f. Soit ℙ une distribution (inconnue) sur 𝒳. L'erreur de m relativement à la cible f est définie par

    ℛ_ℙ,f(m)=ℙ[m(X)≠ f(X)] oùX∼ℙ,

ou écrit de manière équivalente,

    ℛ_ℙ,f(m)=ℙ[{x∈𝒳:m(x)≠ f(x)}].

Pour trouver notre classifieur, il devient nécessaire de supposer qu'il existe un lien entre les données de notre échantillon et le couple (ℙ,f), c'est à dire un modèle de génération des données. On va alors supposer que les x_i sont obtenus par des tirages indépendants suivant ℙ, et qu'ensuite y_i=f(x_i) .

On peut ici définir le risque empirique d'un modèle m,

    ℛ(m)=1/n∑_i=1^n 1(m(x_i)≠ y_i).


Il est alors important d'admettre qu'on ne peut pas trouver un modèle parfait, au sens où ℛ_ℙ,f(m)=0. En effet, si on considère le cas le plus simple qui soit, avec 𝒳={_1,_2} et que ℙ soit telle que ℙ({_1})=p et ℙ({_2})=1-p. Il est aussi possible d'observer x_1 et x_2, et malgré tout, de se tromper sur les étiquettes. Aussi, au lieu de chercher un modèle parfait, on peut tenter d'avoir un modèle approximativement correct. On va alors chercher à trouver m tel que ℛ_ℙ,f(m)≤ϵ, où ϵ est un seuil spécifié a priori.

Une fois admis ce premier écueil, qui fait penser à l'errreur de modèle, notons aussi un second problème. Sur notre exemple à deux valeurs, la probabilité de ne jamais observer x_2 parmi n tirages suivant ℙ est p^n. Il sera alors impossible de trouver m(x_2) car cette valeur n'aura jamais été observée. Autrement dit, aucun algorithme ne peut nous assurer d'avoir avec certitude, avec n observations, ℛ_ℙ,f(m)≤ϵ. On va alors chercher à être probablement approximativement correct (PAC). Pour se faire, on autorise l'algorithme à se tromper avec une probabilité δ, là aussi fixée a priori.

Aussi, quand on construit un classifieur, on ne connaît ni ℙ, ni f, mais on se donne un critère de précision ϵ, et un paramètre de confiance δ, et on dispose de n observations. Notons que n, ϵ et δ peuvent être liés. On cherche alors un modèle m tel que ℛ_ℙ,f(m)≤ϵ avec probabilité (au moins) 1-δ, de manière à être probablement approximativement correct.

Wolpert96 a montré (détaillé dans Wolpert97) qu'il n'existe pas d'algorithme d'apprentissage universel. En particulier, on peut montrer qu'il existe ℙ telle que ℛ_ℙ,f(m) soit relativement grande, avec une probabilité (elle aussi) relativement grande.

L'interprétation qui en est faite est qu'il est nécessaire d'avoir un biais pour apprendre. Comme on ne peut pas apprendre (au sens PAC) sur l'ensemble des fonctions m, on va alors contraindre m à appartenir une classe particulière, notée ℳ. Supposons pour commencer que ℳ contienne un nombre fini de modèles possibles. On peut alors montrer que pour tout ϵ et δ, que pour tout ℙ et f, si on dispose d'assez d'observations (plus précisément n≥ϵ^-1log[δ^-1|ℳ|], alors avec une probabilité plus grande que 1-δ, ℛ_ℙ,f(m^⋆)≤ϵ où

    m^⋆∈m∈ℳargmin{1/n∑_i=1^n 1(m(x_i)≠ y_i)}

autrement dit m^⋆ est un modèle dans ℳ qui minimise le risque empirique.

Un peut aller un peu plus loin, en restant dans le cas où 𝒴={-1,+1}. Une classe ℳ de classifieurs sera dite PAC-apprenable s'il existe n_ℳ:[0,1]^2→ℕ tel que, pour tout ϵ, δ, ℙ et si on suppose que la cible f appartient à ℳ, alors en utilisant n>n_ℳ(ϵ,δ) tirages d'observations x_i suivant ℙ, étiquetés y_i par f, alors il existe m∈ℳ tel que, avec probabilité 1-δ, ℛ_ℙ,f(m)≤ϵ. La fonction n_ℳ est alors appelée complexité d'échantillon pour apprendre. En particulier, nous avons vu que si ℳ contient un nombre fini de classifieurs, alors ℳ est PAC-apprenable avec la complexité n_ℳ(ϵ,δ)=ϵ^-1log[δ^-1|ℳ|].

Naturellement, on souhaiterait avoir un résultat plus général, en particulier si ℳ n'est pas fini. Pour cela, il faut utiliser la dimension VC de Vapnik-Chervonenkis, qui repose sur l'idée de pulvérisation de nuages de points (pour une classification binaire). Considérons k points {x_1,⋯x_k}, et considérons l'ensemble 

    ℰ_k={(m(x_1),⋯,m(x_k)) pour m∈ℳ)}.

Notons que les éléments de ℰ_k appartiennent à {-1,+1}^k. Autrement dit |ℰ_k|≤ 2^k. On dira que ℳ pulvérise l'ensemble des points si toutes les combinaisons sont possibles, c'est à dire |ℰ_k|= 2^k. Intuitivement, les étiquettes de l'ensemble de points ne procurent pas assez d'information sur la cible f, car tout est possible. La dimension VC de ℳ est alors 

    VC(ℳ)=sup{
    k tel que ℳ pulvérise {x_1,⋯x_k}}.



Par exemple si 𝒳=ℝ et que l'on considère l'ensemble des modèles (simples) de la forme m_a,b=1_±(x∈[a,b]). Aucun ensemble de points {x_1,x_2,x_3} ordonnés ne peut être pulvérisé car il suffit d'assigner respectivement +1, -1 et +1 à x_1, x_2 et x_3 respectivement, donc VC<3. En revanche {0,1} est pulvérisé donc VC≥ 2. La dimension de cet ensemble de prédicteur est 2. 
Si on augmente d'une dimension, 𝒳=ℝ^2 et que l'on considère l'ensemble des modèles (simples) de la forme m_a,b=1_±(x∈[a,b]) (où [a,b] désigne le rectangle), alors la dimension de ℳ est ici 4.

Pour introduire les SVM, plaçons nous dans le cas où 𝒳=ℝ^k, et considérons des séparations par des hyperplans passant par l'origine (on dira homogènes), au sens où m_w(x)=1_±(w^ Tx≥ 0). On peut montrer qu'aucun ensemble de k+1 points ne peut être pulvérisé par ces deux espaces homogènes dans ℝ^k, et donc VC(ℳ)=k. Si on rajoute une constante, au sens où m_w,b(x)=1_±(w^ Tx+b≥ 0), on peut montrer qu'aucun ensemble de k+2 points ne peut être pulvérisé par ces deux espaces (non homogènes) dans ℝ^k, et donc VC(ℳ)=k+1.

De cette dimension VC, on déduit le théorème dit fondamental de l'apprentisssage : si ℳ est une classe de dimension d=VC(ℳ), alors il existe des constante positives C et C telles que la complexité d'échantillon pour que ℳ soit PAC-apprenable vérifie

    Cϵ^-1(d+log[δ^-1])
    ≤ n_ℳ(ϵ,δ) ≤Cϵ^-1(dlog[ϵ^-1]+log[δ^-1]).






 §.§ Le choix de l'objectif et la fonction de perte


Ces choix sont essentiels, et dépendent du problème considéré.
Commençons par décrire un modèle historiquement important, le « perceptron » de Rosenblatt, introduit dans des problèmes de classification, où y∈{-1,+1}, inspiré par McCulloch. On dispose de données {(y_i,_i)}, et on va construire de manière itérative un ensemble de modèles m_k(·), où à chaque étape, on va apprendre des erreurs du modèle précédent. Dans le perceptron, on considère un modèle linéaire de telle sorte que :

    m()=1_±(β_0+β≥ 0)={[ +1 si β_0+β≥ 0; -1 si β_0+β< 0 ].,

où les coefficients  sont souvent interprétés comme des « poids » attribués à chacune des variables explicatives. On se donne des poids initiaux (β_0^(0),β^(0)), que l'on va mettre à jour en tenant compte de l'erreur de prédiction commise, entre y_i et la prédiction y_i^(k) :

    y_i^(k)=m^(k)(_i)=1_±(β_0^(k)+β^(k)≥ 0),
 
avec, dans le cas du perceptron :

    β_j^(k+1)=β_j^(k)+η(y-y^(k))^ T_=ℓ(y,y^(k))_j

où ici ℓ(y,y')=1(y ≠ y') est une fonction de perte, qui permettra de donner un prix à une erreur commise, en prédisant y'=m() et en observant y. Pour un problème de régression, on peut considérer une erreur quadratique ℓ_2, telle que ℓ(y,m())=(y - m())^2 ou en valeur absolue ℓ_1, avec ℓ(y,m())=| y - m()|. Ici, pour notre problème de classification, nous utilisions une indicatrice de mauvaise qualification (on pourrait discuter le caractère symétrique de cette fonction de perte, laissant croire qu'un faux positif coûte autant qu'un faux négatif). Une fois spécifiée cette fonction de perte, on reconnaît dans le problème décrit aurapavant une descente de gradient, et on voit que l'on cherche à résoudre :

    m^⋆()=m∈ℳargmin{∑_i=1^n ℓ(y_i,m(_i))
    }

pour un ensemble de prédicteurs ℳ prédéfini. 
Tout problème d'apprentissage machine est mathématiquement formulé comme un problème d'optimisation, dont la solution détermine un ensemble de paramètres de modèle (si la famille ℳ est décrite par un ensemble de paramètres - qui peuvent être des coordonnées dans une base fonctionnelle). On pourra noter ℳ_0 l'espace des hyperplans de ℝ^p au sens où

    m∈ℳ_0    signifie   m()=β_0+   avec   ∈ℝ^p,

engendrant la classe des prédicteurs linéaires. On aura alors l'estimateur qui minimise le risque empirique. Une partie des travaux récents en apprentissage statistique vise à étudier les propriétés de l'estimateur m^⋆, dit « oracle », dans une famille d'estimateurs ℳ,

    m^⋆ = m∈ℳargmin{ℛ(m,m)
    }.

Cet estimateur est, bien entendu, impossible à définir car il dépend de m, le vrai modèle, inconnu.

Mais revenons un peu davantage sur ces fonctions de perte. Une fonction de perte ℓ est une fonction ℝ^d×ℝ^d→ℝ_+, symmétrique, qui vérifie l'inégalité triangulaire, et telle que ℓ(,)=0 si et seulement si =. La norme associée est ·, telle que ℓ(,)=-=ℓ(-,0) (en utilisant le fait que ℓ(,+z)=ℓ(-,z) - nous reverrons cette propriété fondamentale par la suite).

Pour une fonction de perte quadratique, on notera que l'on peut avoir une interprétation particulière de ce problème, puisque :

    y = m∈ℝargmin{∑_i=1^n 1/n [y_i-m]^2 } = m∈ℝargmin{∑_i=1^n  ℓ_2(y_i,m) },

où ℓ_2 est la distance quadratique usuelle
Si l'on suppose - comme on le faisait en économétrie - qu'il existe un modèle probabiliste sous-jacent, et en notant que :

    𝔼(Y) = m∈ℝargmin{‖ Y-m‖^2_ℓ_2} = 
    m∈ℝargmin{𝔼( [Y-m]^2 )  }=
    m∈ℝargmin{𝔼[ℓ_2(Y,m)] }

on notera que ce que l'on essaye d'obtenir ici, en résolvant le problème (<ref>) en prenant pour ℓ la norme ℓ_2, est une approximation (dans un espace fonctionnel donné, ℳ) de l'espérance conditionnelle ↦[Y|=].
Une autre fonction de perte particulièrement intéressante est la perte ℓ_1, ℓ_1(y,m)=| y-m|. Rappelons que

    médiane(y) =  m∈ℝargmin{∑_i=1^n  ℓ_1(y_i,m) }.

Le problème d'optimisation :

    m^⋆=m∈ℳ_0argmin{∑_i=1^n |  y_i-m(_i)|}

est obtenu en économétrie si on suppose que la loi conditionnelle de Y suit une loi de Laplace centrée sur m(), et en maximisant la (log) vraisemblance (la somme des valeurs absolues de erreurs correspond à la log-vraimenblance d'une loi de Laplace). On pourra noter d'ailleurs que si la loi conditionnelle de Y est summétrique par rapport à 0, la médiane et la moyenne conïncident
Si on réécrit cette fonction de perte ℓ_1(y,m)=| (y-m)(1/2-1_y≤ m)|, on peut obtenir une généralisation pour τ∈(0,1) :

    m^⋆_τ=m∈ℳ_0argmin{∑_i=1^n ℓ_τ^ q (y_i,m(_i)) 
    }   avec   ℓ_τ^ q(x,y)= (x-y)(τ-1_x≤ y)

est alors la régression quantile de niveau τ (voir Koenker et dHaultefoeuille). Une autre fonction de perte, introduite par Aigneretal et analysée dans Walltrup, est la fonction associée à la notion d'expectiles :

    ℓ^ e_τ(x,y)= (x-y)^2·|τ-1_x≤ y|

avec τ∈[0,1]. On voit le parallèle avec la fonction quantile :

    ℓ^ q_τ(x,y)= | x-y|·|τ-1_x≤ y|.

KoenkerMachado et YuMoyeed ont d'ailleurs noté un lien entre cette condition et la recherche du maximum de vraisemblance lorsque la loi conditionnelle de Y suit une loi de Laplace assymétrique.

En lien avec cette approche, Gneiting a introduit la notion de « statistique ellicitable » - ou de « mesure  ellicitable » dans sa version probabiliste (ou distributionnelle) : T sera dite « ellicitable » s'il existe une fonction de perte ℓ: ℝ×ℝ→ℝ_+ telle que :

    T(Y)=x∈ℝargmin{∫_ℝℓ(x,y)dF(y)}
    =x∈ℝargmin{𝔼[ ℓ(x,Y)] oùYℒ∼ F
    }.

La moyenne (espérance mathématique) est ainsi ellicitable par la distance quadratique, ℓ_2, alors que la médiane est ellicitable par la distance ℓ_1. Selon Gneiting, cette propriété est essentielle pour construire des prédictions. Il peut alors exister un lien fort entre des mesures associées à des modèles probabilistes et les fonctions de perte. Enfin, la statistique Bayésienne propose un lien direct entre la forme de la loi a priori et la fonction de perte, comme l'ont étudié Berger et Bernardo. Nous reviendrons sur l'utilisation de ces différentes normes dans la section sur la pénalisation.




 §.§ Boosting et apprentissage séquentiel


Nous l'avons vu auparavant: la modélisation repose ici sur la résolution d'un problème d'optimisation, et résoudre le problème décrit par l'équation (<ref>) est d'autant plus complexe que l'espace fonctionnel ℳ est volumineux. L'idée du Boosting,  tel qu'introduit par ShapireFreund, est d'apprendre, lentement, à partir des erreurs du modèle, de manière itérative. À la première étape, on estime un modèle m_1 pour , à partir de , qui donnera une erreur _1. À la seconde étape, on estime un modèle m_2 pour _1, à partir de , qui donnera une erreur _2, etc. On va alors retenir comme modèle, au bout de k itération :

    m^(k)(·)=m_1(·)_∼+m_2(·)_∼_1+m_3(·)_∼_2+⋯+
    m_k(·)_∼_k-1=m^(k-1)(·)+m_k(·).

Ici, l'erreur ε est vue comme la différence entre y et le modèle m(), mais elle peut aussi être vue comme le gradient associé à la fonction de perte quadratique. Formellement,  peut être vu comme un ∇ℓ dans un contexte plus géneral (on retrouve ici une interprétation qui fait penser aux résidus dans les modèles linéaires généralisés).

L'équation (<ref>) peut se voir comme une descente du gradient, mais écrit de manière duale. En effet, la descente de gradient permettant d'obtenir le minimum d'une fonction f repose sur une équation de la forme

    f(x_k)_⟨ f,x_k ⟩∼f(x_k-1)_⟨ f,x_k-1⟩  + (x_k-x_k-1) _α_k∇f(x_k-1)_⟨∇ f,x_k-1⟩

Le problème (<ref>) est dual dans le sens où c'est la fonction f qui doit être optimisée. On pourrait alors écrire une descente de gradient de la forme :

    f_k(x)_⟨ f_k,x⟩∼f_k-1(x)_⟨ f_k-1,x⟩  + (f_k-f_k-1)_β_k⋆_⟨  f_k-1,∇x⟩

où le terme ⋆ peut être interprété comme un gradient, mais dans un espace fonctionnel, et non plus dans ℝ^p. Le problème (<ref>) va alors se réécrire comme un problème d'optimisation :

    m^(k)=m^(k-1)+h∈ℋargmin{∑_i=1^n ℓ(y_i-m^(k-1)(x_i)_ε_k,i,h(x_i))}

où l'astuce consiste à considérer un espace ℋ relativement simple (on parlera de « weak learner » ). Classiquement, les fonctions ℋ sont des fonctions en escalier (que l'on retrouvera dans les arbres de classification et de régression) appelés stumps. Afin de s'assurer que l'apprentissage est effectivement lent, il n'est pas rare d'utiliser un paramètre de « shrinkage », et au lieu de poser, par exemple, ε_1=y-m_1(), on posera ε_1=y-α· m_1() avec α∈[0,1]. On notera que c'est parce qu'on utilise pour ℋ un espace non-linéaire, et que l'apprentissage est lent, que cet algorithme fonctionne bien. Dans le cas du modèle linéaire Gaussien, rappelons en effet que les résidus =- sont orthogonaux aux variables explicatives, , et il est alors impossible d'apprendre de nos erreurs. La principale difficulté est de s'arrêter à temps, car après trop d'itérations, ce n'est plus la fonction m que l'on approxime, mais le bruit. Ce problème est appelé sur-apprentissage.

Cette présentation a l'avantage d'avoir une heuristique faisant penser à un modèle économétrique, en modélisant de manière intérative les résidus par un modèle (très) simple. Mais ce n'est souvent pas la présentation retenue dans la littérature en apprentissage, qui insiste davantage sur une heuristique d'algorithme d'optimisation (et d'approximation du gradient). La fonction est apprise de manière itérative, en partant d'une valeur constante,

    m^(0)=m∈ℝargmin{∑_i=1^n ℓ(y_i,m)
    }.

puis on considère l'apprentissage suivant

    m^(k)=m^(k-1)+h∈ℋargmin∑ _i=1^nℓ(y_i,m^(k-1)(_i)+h(_i)),

qui peut s'écrire, si ℋ est un ensemble de fonctions différentiables,

    m^(k)=m^(k-1)-γ _k∑ _i=1^n∇ _m^(k-1)ℓ(y_i,m^(k-1)(_i)),

où

    γ _k=γargmin ∑ _i=1^nℓ(y_i,m^(k-1)(_i)-γ∇ _m^(k-1)ℓ(y_i,m^(k-1)(_i))).

Pour mieux comprendre le lien avec l'approche décrite auparavant, à l'étape k, on définit des pseudo-résidus en posant

    r_i,k=-.∂ℓ(y_i,m(_i))/∂ m(_i)|_m()=m^(k-1)() pour i=1,⋯,n.

On cherche alors un modèle simple pour expliquer ces pseudo-résidus en fonction des variables explicatives _i, i.e. r_i,k=h^⋆(_i), où h^⋆∈ℋ. Dans un second temps, on cherche un multiplicateur optimal en résolvant

    γ_k = γ∈ℝargmin{∑_i=1^n ℓ(y_i,m^(k-1)(_i)+γ h^⋆(_i))
    }

puis on met à jour le modèle en posant
m_k()=m_k-1()+γ_k h^⋆(). Plus formellement, on passe de l'équation (<ref>) - qui montre clairement qu'on construit un modèle sur les résidus - à l'équation (<ref>) - qui sera ensuite retraduit comme une problème de calcul de gradient - en notant que ℓ(y,m+h)=ℓ(y-m,h). Classiquement, les fonctions ℋ sont construites avec des arbres de régression. Il est aussi possible d'utiliser une forme de pénalisation en posant m_k()=m_k-1()+νγ_k h^⋆(), avec ν∈(0,1). Mais revenons un peu plus longuement sur l'importance de la pénalisation avant de discuter les aspects numériques de l'optimisation.



 §.§ Pénalisation et choix de variables


Dans la section <ref>, nous avions évoqué le principe de parcimonie, populaire en économétrie. Le critère d'Akaike était basé sur une pénalisation de la vraisemblance en tenant compte de la complexité du modèle (le nombre de variables explicatives retenues). Si en économétrie, il est d'usage de maximiser la vraisemblance (pour construire un estimateur asymptotiquement sans biais), et de juger de la qualité du modèle ex-post en pénalisant la vraisemblance, la stratégie ici sera de pénaliser ex-ante dans la fonction objectif, quitte à construire un estimateur biaisé. Typiquement, on va construire :

    (β_0,λ,_λ)=argmin{∑_i=1^n ℓ(y_i,β_0+)+λ pénalisation( β)
    },

où la fonction de pénalisation sera souvent une norme ‖·‖ choisie a priori, et un paramètre de pénalisation λ (on retrouve en quelque sorte la distinction entre AIC et BIC si la fonction de pénalisation est la complexité du modèle - le nombre de variables explicatives retenues). Dans le cas de la norme ℓ_2, on retrouve l'estimateur Ridge, et pour la norme ℓ_1, on retrouve l'estimateur lasso (« Least Absolute Shrinkage and Selection Operator ⁠»). La pénalisation utilisée auparavant faisait intervenir le nombre de degrés de liberté du modèle, il peut alors paraître surprenant de faire intervenir ‖β‖_ℓ_2 comme dans la régression Ridge. On peut toutefois envisager une vision Bayésienne de cette pénalisation. Rappelons que dans un modèle Bayésien :

    ℙ[θ|y]_a posteriori∝ℙ[y|θ]_vraisemblance·ℙ[θ]_a priori   soit   logℙ[θ|y]= logℙ[y|θ]_log vraisemblance + logℙ[θ]_pénalisation.

Dans un modèle linéaire Gaussien, si on suppose que la loi a priori de θ suit une loi normale centrée, on retrouve une pénalisation basée sur une forme quadratique des composantes de θ.



Avant de revenir en détails sur ces deux estimateurs, obtenus en utilisant la norme ℓ_1 ou la norme ℓ_2, revenons un instant sur un problème très proche : celui du meilleur choix de variables explicatives. Classiquement (et ça sera encore plus vrai en  grande dimension), on peut disposer d'un grand nombre de variables explicatives, p, mais beaucoup sont juste du bruit, au sens où β_j=0 pour un grand nombre de j. Soit s le nombre de covariables (réellement) pertinentes, s=#𝒮
avec 𝒮()={ j=1,⋯,p; β_j≠ 0}. Si on note _𝒮 la matrice constituée des variables pertinentes (en colonnes), alors on suppose que le vrai modèle
est de la forme y=_𝒮_𝒮+ε. Intuitivement, un estimateur intéressant serait alors β_𝒮=[_𝒮_𝒮]^-1_𝒮, mais cet estimateur n'est que théorique car 𝒮 est ici inconnue. Cet estimateur est l'estimateur oracle évoqué auparavant.
On peut alors être tenté de résoudre 

    (β_0,s,_s)=argmin{∑_i=1^n ℓ(y_i,β_0+)}, sous la contrainte #𝒮()=s.


Ce problème a été introduit par FosterGeorge94 en introduisant la norme ℓ_0. Plus précisément, définissons ici les trois normes suivantes

    ‖a‖_ℓ_0=
    ∑_i=1^d 1(a_i≠ 0),   ‖a‖_ℓ_1=
    ∑_i=1^d |a_i|   et   ‖a‖_ℓ_2=(
    ∑_i=1^d a_i^2)^1/2, pour a∈ℝ^d.








Considérons les problèmes d'optimisation de la Table <ref>. Si on considère le problème classique où ℓ est la norme quadratique, les deux problèmes de l'équation (ℓ1) de la Table <ref> sont équivalents, au sens où, pour toute solution (β^⋆,s^⋆) au problème de gauche, il existe λ^⋆ tel que (β^⋆,λ^⋆) soit solution du problème de droite; et inversement. Le résultat est également vrai pour les problèmes (ℓ2)[ Pour (ℓ1), s'il y a équivalence au niveau théorique, il peut exister des soucis numériques car il n'y a pas forcément unicité de la solution.]. Il s'agit en effet de problèmes convexes. En revanche, les deux problèmes (ℓ0) ne sont pas équivalents : si pour (β^⋆,λ^⋆) solution du problème de droite, il existe s^⋆ tel que β^⋆ soit solution du problème de gauche, la récriproque n'est pas vraie. Plus généralement, si on veut utiliser une norme ℓ_p, la sparsitée est obtenue si p≤ 1 alors qu'il faut avoir p≥ 1 pour avoir la convexité du programme d'optimisation.

On peut être tenté de résoudre le programme pénalisé (ℓ0) directement, comme le suggère FosterGeorge94. Numériquement, c'est un problème combinatoire complexe en grande dimension (Natarajan note que c'est un problème NP-difficile), mais il est possible de montrer que si λ∼σ^2 log(p), alors

    𝔼([-_0]^2) ≤𝔼(_𝒮_𝒮-_0]^2)_=σ^2 #𝒮·(4log p+2+o(1)).

Notons quand dans ce cas

    _λ,j^ sub = {[       0  si  j∉𝒮_λ(); _j^ ols si  j∈𝒮_λ(), ].

où 𝒮_λ() désigne l'ensemble des coordonnées non nulle lors de la résolution de (ℓ0).

Le problème  (ℓ2) est strictement convexe si ℓ est la norme quadratique, autrement dit, l'estimateur Ridge est toujours bien défini, avec en plus une forme explicite pour l'estimateur, 

    β_λ^ ridge=(+λ𝕀)^-1
    =(+λ𝕀)^-1()β^ ols.

Aussi, on peut en déduire que

    biais[β_λ^ ridge]=-λ[+λ𝕀]^-1 β^ ols et Var[β_λ^ ridge]=
    σ^2[+λ𝕀]^-1[+λ𝕀]^-1.

Avec une matrice de variables explicatives orthonormées (i.e. =𝕀), les expressions se simplifient

    biais[β_λ^ ridge]=λ/1+λ β^ ols et Var[β_λ^ ridge]=σ^2/(1+λ)^2𝕀.

Notons que Var[β_λ^ ridge]<Var[β^ ols]. En notant que

    mse[β_λ^ ridge]=kσ^2/(1+λ)^2+λ^2/(1+λ)^2,

on obtient une valeur optimale pour λ: λ^⋆=kσ^2 /.


En revanche, si ℓ n'est plus la norme quadratique mais la norme ℓ_1, le problème  (ℓ1) n'est pas toujours strictement convexe, et en particulier, l'optimum n'est pas toujours unique (par exemple si  est singulière). Mais le fait que ℓ soit strictement convexe β sera unique. Notons de plus que deux solutions sont forcément cohérentes en terme de signe des coefficients : il n'est pas possible d'avoir β_j<0 pour une solution et β_j>0 pour une autre. D'un point de vue heuristique, le programme (ℓ1) est intéressant car il permet d'obtenir dans bon nombre de cas une solution en coin, qui correspond à une résolution de problème de type (ℓ0) - comme le montre de manière visuelle la Figure <ref>. 



Considérons un modèle très simple: y_i=x_iβ+ε, avec une pénalité ℓ_1 et une fonction de perte ℓ_2. Le problème (ℓ2) sécrit alors

    min{-2β+ββ+2λ|β|}

La condition du premier ordre est alors

    -2 + 2β± 2λ=0.

le signe du dernier terme dépend du signe de β.
Supposons que l'estimateur par moindre carrés (obtenu en posant λ=0) soit (strictement positif), autrement dit >0. Si λ n'est pas trop grand, on peut imaginer que β soit du même signe que β^ mco, et donc la condition devient

    -2 + 2β+ 2λ=0.

et la solution est 

    β_λ^ lasso=-λ/.

En augmentant λ, on va arriver à un moment où β_λ=0. Si on augmente encore un peu β_λ ne devient pas négatif car dans ce cas le dernier terme de la condition du premier ordre change, et dans ce cas on cherche à résoudre

    -2 + 2β- 2λ=0.

dont la solution est alors

    β_λ^ lasso=+λ/.

Mais cette solution est positive (nous avions supposé >0), et donc il est possible d'avoir en même temps β_λ<0. Aussi, au bout d'un moment, β_λ=0, qui est alors une solution de coin. Les choses sont bien entendu plus compliquées en dimension plus grande (TibWasserman revient longuement sur la géométrie des solutions) mais comme le note Candes, sous des hypothèses minimales garantissant que les prédicteurs ne sont pas fortement corrélées, le Lasso obtient une erreur quadratique presque aussi bonne que si l'on dispose d'un oracle fournissant des informations parfaites sur quels β_j sont non nulles.
Moyennant quelques hypothèses techniques supplémentaires, on peut montrer que cet estimateur est « sparsistent »  au sens où le support de β_λ^ lasso est celui de , autement dit Lasso a permis de faire de la sélection de variables (plus de discussions sur ce point peuvent être obtenues dans HTT).


De manière plus générale, on peut montrer que β_λ^ lasso est un estimateur biaisé, mais qui peut être de variance suffisamment faible pour que l'erreur quadratique moyenne soit plus faible qu'en utilisant des moindres carrés. Pour comparer les trois techniques, par rapport à l'estimateur par moindre carrées (obtenu quand λ=0), si on suppose que les variables explicatives sont orthonormées, alors

    β_λ,j^ sub=β_j^ ols1_|β_λ,j^ sub|>b,   β_λ,j^ ridge=β_j^ ols/1+λ  et  β_λ,j^ lasso=signe[β_j^ ols]·(|β_j^ ols|-λ)_+.






 §.§ Optimisation et aspects algorithmiques


En économétrie, l'optimisation (numérique) est devenu omniprésente dès que l'on a quitté le modèle Gaussien. Nous l'avions rapidement évoqué dans la section sur la famille exponentielle, et l'utilisation du score de Fisher (descente de gradient) pour résoudre la condition du premier ordre X^TW(β)^-1[y-]=0. En apprentissage, l'optimisation est l'outil central. Et il est nécessaire d'avoir des algorithmes d'optimisation efficaces, pour résoudre des problèmes de la forme :

    ∈∈ℝ^pargmin{∑_i=1^n ℓ(y_i,β_0+)+λ‖β‖}.


Dans certains cas, au lieu de faire de l'optimisation globale, il suffit de considérer de l'optimisation par coordonnées (largement étudiée dans Daubechies). Si f:ℝ^d→ℝ est convexe et différentiable, alors 

    si  vérifie f(+he_i)≥ f() pour tout h>0 et i∈{1,⋯, d}, alors 
    f()=min{f},

où e=(e_i) est la base canonique de  ℝ^d. Cette propriété n'est toutefois pas vraie dans le cas non-différentiable. Mais si on suppose que la partie non-différentiable est séparable (additivement), elle redevient vraie. Plus précisément, si 

    f()=g()+∑_i=1^d h_i(x_i) avec {[ g: ℝ^d→ℝ convexe-différentiable;               h_i: ℝ→ℝ convexe. ].

C'est le cas pour la régression Lasso, f()=-X_ℓ_2+λ_ℓ_1, comme le montre Tsen.
On peut alors utiliser un algorithme de descente par coordonnées: à partir d'une valeur initiale ^(0), on considère (en itérant)

    x_j^(k)∈argmin{
    f(x_1^(k),⋯,x_k-1^(k),x_k,x_k+1^(k-1),
    ⋯,x_n^(k-1))
    }  pour j=1,2,⋯,n.

Ces problèmes algorithmiques peuvent paraître secondaires à des économètres. Ils sont pourtant essentiels en apprentissage machine: une technique est intéressante s'il existe un algorithme stable et rapide, qui permet d'obtenir une solution. 
 Ces techniques d'optimisation sont d'ailleurs transposables : par exemple, on pourra utiliser cette technique de descente par coordonnées dans le cas des méthodes svm (dit « à support vecteur⁠») lorsque l'espace n'est pas linéairement séparable, et qu'il convient de pénaliser l'erreur de classification (nous reviendrons sur cette technique dans la prochaine section).




 §.§ In-sample, out-of-sample et validation croisée


Ces techniques semblent intellectuellement intéressantes, mais nous n'avons pas encore abordé le choix du paramètre de pénalisation λ. Mais ce problème est en fait plus général, car comparer deux paramètres _λ_1 et _λ_2 revient en fait à comparer deux modèles. En particulier, si on utilise une méthode de type Lasso, avec des seuils λ différents, on compare des modèles qui n'ont pas la même dimension. Dans la section <ref>, nous avions abordé le problème de la comparaison de modèles sous l'angle économétrique (en pénalisant les modèles trop complexes). Dans la littérature en apprentissage, juger de la qualité d'un modèle sur les données qui ont servi à le construire ne permet en rien de savoir comment le modèle se comportera sur des nouvelles données. Il s'agit du problème dit de « généralisation ». L'approche classique consiste alors à séparer l'échantillon (de taille n) en deux : une partie qui servira à entraîner le modèle (la base d'apprentissage, in-sample, de taille m) et une partie qui servira à tester le modèle (la base de test, out-of-sample, de taille n-m). Cette dernière permet alors de mesure un vrai risque prédictif. Supposons que les données soient générées par un modèle linéaire y_i=_i_0 + ε_i où les ε_i sont des réalisations de lois indépendantes et centrées. Le risque quadratique empirique in-sample est ici

    1/m∑_i=1^m𝔼([_i-_i_0]^2)=𝔼([_i-_i_0]^2),

pour n'importe quelle observation i. En supposant les résidus ε Gaussiens, alors on peut montrer que ce risque vaut σ^2trace(Π_𝒳)/m soit σ^2p/m. En revanche le risque quadratique empirique out-of-sample est ici

    𝔼([-_0]^2)

où  est une nouvelle observation, indépendante des autres. On peut noter que 

    𝔼([-_0]^2|)=
    
    σ^2()^-1,

et en intégrant par rapport à ,

    𝔼([-_0]^2)=
    𝔼(𝔼([-_0]^2|))=
    σ^2trace(𝔼[]𝔼[()^-1]).

L'expression est alors différente de celle obtenue in-sample, et en utilisation la majoration de GrovesRothenberg69, on peut montrer que

    𝔼([-_0]^2) ≥σ^2p/m,

ce qui est assez intuitif, finalement. Hormis certains cas simple, il n'y a pas de formule simple. Notons toutefois que si ∼𝒩(0,σ^2𝕀), alors  suit une loi de Wishart, et on peut montrer que

    𝔼([-_0]^2)=σ^2p/m-p-1.

Si on regarde maintenant la version empirique: si β est estimé sur les m premières observations,

    ℛ^  IS=∑_i=1^m [y_i-x_iβ]^2
       et   ℛ^ OS=∑_i=m+1^n [y_i-x_iβ]^2,

et comme l'a noté Leeb, ℛ^  IS -ℛ^  OS≈ 2 ·ν où ν représente le nombre de degrés de libertés, qui n'est pas sans rappeler la pénalisation utilisée dans le critère d'Akaike.

La Figure <ref> montre l'évolution respective de ℛ^ IS et ℛ^ OS en fonction de la complexité du modèle (nombre de degrés dans une régression polynomiale, nombre de noeuds dans des splines, etc). Plus le modèle est complexe, plus ℛ^ IS va diminuer (c'est la courbe rouge). Mais ce n'est pas ce qui nous intéresse ici : on veut un modèle qui prédise bien sur de nouvelles données (autrement dit out-of-sample). Comme le montre la Figure <ref>, si le modèle est trop simple, il prédit mal (tout comme sur les données in-sample). Mais ce que l'on peut voir, c'est que si le modèle est trop complexe, on est dans une situation de « sur-apprentissage » : le modèle va commencer à modéliser le bruit.




Au lieu de séparer la base en deux, avec une partie des données qui vont servir à calibrer le modèle et une autre à étudier sa performance, il est aussi possible d'utiliser la validation croisée. Pour présenter l'idée générale, on peut revenir au « jackknife », introduit par Quenouille1 (et formalisé par Quenouille2 et Tukey) et utilisé en statistique pour réduire le biais. En effet, si on suppose que { y_1,⋯,y_n} est un échantillon tiré suivant une loi F_θ, et que l'on dispose d'un estimateur T_n()=T_n(y_1,⋯,y_n), mais que cet estimateur est biaisé, avec
[T_n()]=θ+O(n^-1),  
il est possible de réduire le biais en considérant :

    T_n()=1/n∑_i=1^n T_n-1(_(i))   avec   _(i)=(y_1,⋯,y_i-1,y_i+1,⋯,y_n).

On peut alors montrer que [T_n()]=θ+O(n^-2).

L'idée de la validation croisée repose sur l'idée de construire un estimateur en enlevant une observation. Comme on souhaite construire un modèle prédictif, on va comparer la prévision obtenue avec le modèle estimé, et l'observation manquante :

    ℛ^ CV=1/n∑_i=1^n ℓ(y_i,m_(i)(_i))

On parlera ici de méthode « leave-one-out » (loocv).

On utilise classiquement cette technique pour trouver le paramètre optimal dans les méthodes de lissage exponentiel, pour des séries chronologiques. Dans le lissage simple, on va construire une prédiction de la forme _ty_t+1=α·_t-1y_t+(1-α)· y_t, avec α∈[0,1], et on va considérer :

    α^⋆ = α∈[0,1]argmin{∑_t=2^T ℓ(_t-1y_t,y_t) },

comme le décrit Hyndman.

Le principal problème de la méthode « leave-one-out » est qu'elle nécessite de calibrer n modèles, ce qui peut être problématique en grande dimension. Une méthode alternative est la validation croisée par k-blocs (dit « k-fold cross validation ») qui consiste à utiliser une partition de {1,⋯,n} en k groupes (ou blocs) de même taille,
ℐ_1,⋯,ℐ_k, et notons ℐ_j={
1,⋯,n}\ℐ_j. En notant m_(j) construit sur l'échantillon ℐ_j, on pose alors :

    ℛ^k- CV=1/k∑_j=1^k ℛ_j   où  ℛ_j=k/n∑_i∈ℐ_jℓ(y_i,m_(j)(_i)).

La validation croisée standard, où une seule observation est enlevée à chaque fois (loocv), est un cas particulier, avec  k=n. Utiliser k=5,10 a un double avantage par rapport à k=n : (1) le nombre d'estimations à effectuer est beaucoup plus faible, 5 ou 10 plutôt que n ; (2) les échantillons utilisés pour l'estimation sont moins similaires et donc, moins corrélés les uns aux autres, ce qui tend à éviter les excès de variance, commme le rappelle James.











Une autre alternative consiste à utiliser des échantillons boostrappés. Soit ℐ_b un échantillon de taille n obtenu en tirant avec remise dans {1,⋯,n} pour savoir quelles observations (y_i,_i) seront gardées dans la population d'apprentissage (à chaque tirage). Notons ℐ_b={
1,⋯,n}\ℐ_b. En notant m_(b) construit sur l'échantillon ℐ_b, on pose alors :

    ℛ^ B=1/B∑_b=1^B ℛ_b   où  ℛ_b=n_b/n∑_i∈ℐ_bℓ(y_i,m_(b)(_i)),

où n_b est le nombre d'observations qui n'ont pas été conservées dans ℐ_b. On notera qu'avec cette technique, en moyenne e^-1∼ 36.7% des observations ne figurent pas dans l'échantillon boostrappé, et on retrouve un ordre de grandeur des proportions utilisées en créant un échantillon de calibration, et un échantillon de test.
 En fait, comme l'avait montré Stone, la minimization du AIC est à rapprocher du critère de validation croisée, et Shao97 a montré que la minimisation du BIC correspond à de la validation croisée de type k-fold, avec k=n/log n.



§ QUELQUES OUTILS DE MACHINE LEARNING





 §.§ Réseaux de Neurones




Les réseaux de neurones sont des modèles semi-paramétriques. Néanmoins, cette famille de modèles peut être appréhendée de la même manière que les modèles non-paramétriques: la structure des réseaux de neurones (présentée par la suite) peut être modifiée afin d'étendre la classe des fonctions utilisées pour approcher une variable d'intérêt. Plus précisément, Cybenko a démontré que l'ensemble des fonctions neuronales est dense dans l'espace des fonctions continues sur un compact. En d'autres termes, on a un cadre théorique permettant de garantir une forme d'approximation universelle. Il impose en outre une définition d'un neurone et met en avant l'existence d'un nombre de neurones suffisant pour approcher toute fonction continue sur un compact. Ainsi, un phénomène continue peut être approché par une suite de neurones: on appellera cette suite « réseau de neurones à une couche ». Si ce théorème d'approximation universelle est démontré en 1989, le premier neurone artificiel fonctionnel fut introduit par Franck Rosenblatt au milieu du XXième siècle, dans Rosenblatt. Ce neurone, qualifié de nos jours de « neurone élémentaire », porte le nom de « Perceptron ». Il a permis dans ses premières utilisations de déterminer le sexe d'un individu présenté aux travers d'une photo. Si ce premier neurone est important, c'est qu'il introduit le premier formalisme mathématique d'un neurone biologique. On peut décrire un neurone artificiel par analogie avec une cellule nerveuse : 



  -  les synapses apportant l'information à la cellule sont formalisés par un vecteur réel. La dimension du vecteur d'entrée du neurone (qui n'est d'autre qu'une fonction) correspond biologiquement au nombre de connections synaptiques;

  -  chaque signal apporté par un synapse est ensuite analysé par la cellule. Mathématiquement, ce schéma est transcrit par la pondération des différents éléments constitutifs du vecteur d'entrée;

  - en fonction de l'information acquise, le neurone décide de retransmettre ou non un signal. Ce phénomène est répliqué par l'introduction d'une fonction d'activation. Le signal de sortie est modélisé par un nombre réel calculé comme image par la fonction d'activation du vecteur d'entrée pondéré.



Ainsi, un neurone artificiel est un modèle semi-paramétrique. Le choix de la fonction d'activation est en effet laissé à l'utilisateur. Nous introduisons dans le paragraphe qui suit une formalisation rigoureuse qui nous permettra de poser le modèle, et de faire le lien avec les notations économétriques usuelles. On peut alors définir un neurone élémentaire formellement par : 




  -  un espace d'entrée 𝒳, généralement ℝ^k avec k ∈ℕ^*;

  -  un espace de sortie 𝒴, généralement ℝ ou un ensemble fini (classiquement {0,1}, mais on préférera ici {-1,+1});

  -  un vecteur de paramètres w∈ℝ^p 

  -  une fonction d'activation ϕ : ℝ→ℝ.  Cette fonction doit être dans l'idéal monotone, dérivable et bornée (on dira ici « saturante  ») afin de s'assurer de certaines propriétés de convergence. 

Cette dernière fonction ϕ fait penser aux transformations logistique ou probit, populaire en économétrie (qui sont des fonctions de répartition, à valeur dans [0,1], idéal quand 𝒴 est l'ensemble {0,1}). Pour les réseaux de neurones, on utilisera plutôt la tangente hyperbolique, la fonction arctangente ou les fonctions sigmoïdes pour des problèmes de classification (𝒴={-1,+1}.
On appellera neurone toute application f_w de 𝒳 dans 𝒴 définie par : 

    y=f_w() = ϕ (wx),   ∀x∈𝒳.
















Pour le perceptron introduit par Rosenblatt, on assimile un neurone élémentaire à la fonction :

    y=f_w(x)=signe(wx)  ∀x∈𝒳
 

On remarque que selon cette formalisation, beaucoup de modèles statistiques comme par exemple les régressions logistiques pourraient être vus comme des neurones. En effet si l'on regarde d'un peu plus près, tout modèle glm (« Generalized Linear Model ») pourrait s'interpréter comme un neurone formel où la fonction d'activation ϕ n'est d'autre que l'inverse de la fonction de lien canonique (par exemple). Si g désigne la fonction de lien du glm , w le vecteur de paramètres, y la variable à expliquer et x le vecteur des variables explicatives de même dimension que w : 

    g(𝔼[Y|X=x])=wx
 
On retrouve la modélisation neuronale en prenant ϕ = g^-1. Cependant, là où réside la différence majeure entre les glm et le modèle neuronale est que ce dernier n'introduit aucune hypothèse de distribution sur Y|X (on n'a d'ailleurs pas besoin d'introduire ici de modèle probabiliste). D'autre part, lorsque le nombre de neurones par couche augmente, la convergence n'est pas nécessairement garantie si la fonction d'activation ne vérifie pas certaines propriétés (qu'on ne retrouve pas dans la majorité des fonctions de liens canoniques des glm). 
Cependant, comme énoncé précédemment, la théorie des réseaux de neurones introduit des contraintes mathématiques supplémentaires sur la fonction g (détaillé dans Cybenko). Ainsi par exemple, une régression logistique peut être perçue comme un neurone alors que les régressions linéraires généralisées ne vérifient pas toutes les hypothèses nécessaires.



Toujours par analogie avec le fonctionnement du système nerveux, il est alors possible de connecter différents neurones entre eux. On parlera de structure de réseaux de neurones par couche. Chaque couche de neurones recevant à chaque fois le même vecteur d'observation. 

Pour revenir à une analogie plus économétrique, on peut imaginer passer par une étape intermédiaire (on reviendra sur cette construction dans la Figure <ref>), par exemple en ne faisant pas une régression sur les variables brutes  mais un ensemble plus faible de variables orthogonales, obtenues par exemple suite à une analyse en composantes principales. Soit A la matrice associée à cette transformation linéaire, avec A de taille k× p si on souhaite utiliser les p premières composantes. Notons z la transformation de x, au sens où z = A, ou encore z_j=A_j. Un généralisation du modèle précédant peut être de poser

    y=f(x) = ϕ (wz)= ϕ (wA)=,   ∀x∈𝒳,

où cette fois w∈ℝ^p. On a ici une transformation linéaire (en considérant une analyse en composante principale) mais on peut imaginer une généralisation avec des transformée non-linéaire, avec une fonction de la forme

    y=f(x) = ϕ (w F_A()=,   ∀x∈𝒳,

où F est ici une fonction ℝ^k→ℝ^p. C'est le réseau de neurone à deux couches. Plus généralement, pour formaliser la construction, on introduit les notations suivantes : 



  - K ∈ℕ^* : nombre de couches;

  - ∀ k ∈{ 1 ,⋯ K } , p_k représente le nombre de neurones dans la couche k;

  - ∀ k ∈{ 1 ,⋯ K } ,W_k désigne la matrice des paramètres associés à la couche k. Plus précisément, W_k est une matrice p_k× p_k-1 et pour tout l ∈{ 1 ,⋯ p_k } , w_k,l∈ℝ^ p_k-1 désigne le vecteur de poids associé au neurone élémentaire l de la couche k;

  - on appellera W = {W_1,..,W_K}, l'ensemble des paramètres associés au réseau de neurones.

  - F^k_W_k : ℝ^p_k-1→ℝ^p_k désigne la fonction de transfert associé à la couche k. Pour des raisons de simplification, on pourra également écrire F^k;

  - ŷ_k ∈ℝ^p_k représentera le vecteur image de la couche k ∈{ 1 ,⋯, K };

  - on appellera F= F_W = F^1 ∘⋯∘ F^K la fonction de transfert associée au réseau global. A ce titre, si x∈𝒳, on pourra noter y= F_W(x).





La Figure <ref> permet d'illustrer les notations présentées ici[Source: http://intelligenceartificielle.org.]. Chaque cercle représente un neurone élémentaire. Chaque rectangle englobant plusieurs cercles représente une couche. On parle de couche d'entrée pour la première couche prenant en « input » les observation x∈𝒳, de couche de sortie pour la couche fournissant en « output » la prédiction ŷ∈𝒴. Les autres couches sont couramment appelées couches cachées. 

Un réseau de neurones multicouches est donc également un modèle semi-paramétrique dont les paramètres sont l'ensemble des composantes des matrices W_k pour tout entier k de { 1,⋯, K }. Chaque fonction d'activation associée à chaque neurone (chaque cercle de la Figure <ref>) est à déterminer par l'utilisateur. 




Une fois que les paramètres à calibrer du modèle sont identifiés (ici les réels constituant les matrices W_k pour chaque couche k ∈{1,⋯,K }) , il est nécessaire de fixer une fonction de perte ℓ. En effet, on rappelle que l'objectif de l'apprentissage supervisé sur une base d'apprentissage de n ∈ℕ^* couples (y_i,x_i) ∈𝒴×𝒳 est de minimiser le risque empirique : 

    ℛ_n(F_W) = 1/n∑_i=1^n ℓ(y_i,F_W(x_i))

        
Afin d'illustrer les propos, intéressons nous à l'exemple suivant qui illustrera également la démarche opérée. Supposons que nous observons un phénomène y aux travers de n observations y_i ∈ [-1,1]. On souhaiterait expliquer ce phénomène à partir des variables explicatives  que l'on suppose à valeurs réelles. La « théorie de l'approximation universelle » nous indique qu'un réseau à une couche de neurones devrait permettre de modéliser le phénomène (sous hypothèse qu'il soit continue). On note toutefois que ce théorème ne donne pas de vitesse de convergence. Il est alors laissé à l'utilisateur le choix de la structure. Ainsi par exemple une première structure pourrait être un simple neurone dont la fonction d'activation serait la fonction tangente hyperbolique.

On aurait ainsi comme premier modèle:

    y_1 = tanh(w_0 + w_1 x)
 
où les paramètres w_0, w_1  sont les paramètres à optimiser de sorte que sur les données d'apprentissage, le risque empirique soit minimal. 

 Si l'on suit toujours la philosophie du théorème d'approximation universelle, en ajoutant plusieurs neurones, l'erreur est censée diminuer. Cependant,  ne connaissant pas la fonction à estimer, on ne peut l'observer qu'aux travers de l'échantillon. Ainsi, mécaniquement, on s'attend  à ce que plus on ajoute de paramètres, plus l'erreur sur la base d'apprentissage diminue. L'analyse de l'erreur sur la base de test permet alors d'évaluer notre capacité à généraliser (cf partie précédente).
 
On peut ainsi s'intéresser à un second modèle qui cette fois utilise plusieurs neurones. Par exemple, considérons le modèle

    y_2 = w_a tanh(w_0 + w_1 x) + w_b tanh(w_2 + w_3 x) + w_c tanh(w_4 + w_5 x)
 
où les paramètres w_0,..,w_5 ainsi que w_a, w_b, w_c  sont les paramètres à optimiser. 
Calibrer un réseaux de neurones revient alors à réitérer ces étapes de modification de la structure jusqu'à minimisation du risque sur la base de test. 
 









Pour une structure de réseau de neurones fixée (c'est-à-dire nombre de couches, nombre de neurones par couches et fonctions d'activation fixés), le programme revient donc à déterminer l'ensemble de paramètres W^*=(W_1,...,W_K) de sorte que :

    W^*∈ W=(W_1,...,W_K)argmin{1/n∑_i=1^n ℓ(y_i,F_W(_i))
    }.

De cette formule apparaît l'importance du choix de la fonction ℓ. Cette fonction de perte, elle quantifie l'erreur moyenne commise par notre modèle F_W sur la base d'apprentissage. A priori ℓ peut être choisie arbitrairement. Cependant, dans l'optique de résoudre un programme d'optimisation, on préféra des fonctions de coût sous-différentiables et convexes afin de garantir la convergence des algorithmes d'optimisation. 
Parmi les fonctions de perte classiques, en plus de la fonction de perte quadratique ℓ_2 on retiendra la fonction dite « Hinge » - ℓ(y,ŷ)=max(0, 1-yŷ) - ou la fonction dite logistique - ℓ(y,ŷ)=log(1-e^-yŷ).

















En définitive les réseaux de neurones sont des modèles semi-paramétriques dont le nombre de paramètres est croissant avec le nombre de couches et de neurones par couche. Il est laissé à l'utilisateur de choisir les fonctions d'activation et la structure du réseau. Ceci explique l'analogie avec la philosophie des modèles non-paramétriques faite auparavant.



Les réseaux de neurones ont été utilisés très tôt en économie et en finance, en particulier sur les défauts d'entreprises - Tam ou Altman - ou plus récemment la notation de crédit - Blanco ou Khashman. Cependant les structures telles que présentées précédemment sont généralement limitées. L'apprentissage profond (ou « deep learning ») caractérise plus particulièrement des réseaux de neurones plus complexes (parfois plus d'une dizaine de couches avec parfois des centaines de neurones par couche). Si aujourd'hui ces structures sont très populaires en analyse du signal (image, texte, son) c'est qu'elles sont capables à partir d'une quantité d'observations très importante d'extraire des informations que l'humain ne peut percevoir et de faire face à des problèmes non linéaires, comme le rappelle LeCun. 

L'extraction d'informations peut, par exemple, se faire grâce à la convolution. Procédé non supervisé, il a permis notamment d'obtenir d'excellente performance dans l'analyse d'image. Techniquement, cela peut s'apparenter à une transformation à noyaux (comme utilisé dans les techniques SVM). Si une image peut être perçue comme une matrice dont chaque coordonnée représente un pixel, une convolution reviendrait à appliquer une transformation sur un point (ou une zone) de cette matrice générant ainsi une nouvelle donnée. Le procédé peut ainsi être répété en appliquant des transformations différentes (d'où la notion de couches convolutives). Le vecteur final obtenu peut alors enfin alimenter un modèle neuronal comme introduit dans le paragraphe précédant. En fait, plus généralement, une couche de convolution peut être perçue comme un filtre qui permet de transformer la donnée initiale.










Une explication intuitive pour laquelle l'apprentissage approfondi, en particulier les réseaux nerveux profonds, est si puissant pour décrire des relations complexes dans les données, c'est leur construction autour de l'approximation fonctionnelle simple et l'exploitation d'une forme de hiérarchie, comme le note Lin.
Néanmoins les modèles de type « deep learning » sont plus difficiles à appréhender car ils nécessitent beaucoup de jugement empirique. En effet, si aujourd'hui les bibliothèques open sources (keras, torch, etc.) permettent de paralléliser plus facilement les calculs en utilisant par exemple les GPU (Graphical Processor Units), il reste néanmoins à l'utilisateur de déterminer la structure du réseau de neurones le plus approprié.




 §.§ Support Vecteurs Machine


Comme nous l'avions noté auparavant, dans les problèmes de classification en apprentissage machine (comme en traitement du signal) on préférera avoir des observations dans l'ensemble {-1,+1} (plutôt que {0,1}, comme en économétrie). Avec cette notation, Cortes ont posé les bases théorique des modèles dit svm, proposant une alternative aux réseaux de neurones alors très populaires comme algorithme de classification dans la communauté de l'apprentissage machine. L'idée initiale des méthodes de « Support Vectors Machine » (svm) consiste à trouver un hyperplan séparateur divisant l'espace en deux ensembles de points le plus homogène possible (i.e. contenant des labels identiques). En dimension deux, l'algorithme consiste à déterminer une droite séparant l'espace en deux zones les plus homogènes possibles. La résolution de ce problème possédant parfois une infinité de solution (il peut en effet exister une infinité de droites qui séparent l'espace en deux zones distinctes et homogènes), on rajoute généralement une contrainte supplémentaire. L'hyperplan séparateur doit se trouver le plus éloigné possible des deux sous-ensembles homogènes qu'il engendre. On parlera ainsi de marge. L'algorithme ainsi décrit est alors un svm linéaire à marge.



Si un plan peut être caractérisé entièrement par un vecteur directeur w orthogonal à ce dernier et une constante b, appliquer un algorithme SVM à un ensemble de n∈ℕ^* points _i de ℝ^p labellisés par y_i ∈{-1,1} revient alors à résoudre un programme d'optimisation sous contrainte similaire à celui d'un lasso (distance quadratique sous contrainte linéaire). Plus particulièrement, on sera amené à résoudre :

    (w^⋆,b^⋆) = w,bargmin{w^2}=w,bargmin{ww},


    sous contrainte  ∀ i ∈{ 1, ⋯, n }, {[ ω^ Tx_i+b ≥ +1   lorsque   y_i=+1; ω^ Tx_i+b ≤ -1   lorsque   y_i=-1;                                   ].
 

La contrainte peut être relâchée en autorisant que dans un sous-ensemble, un point puisse ne pas être du même label que la majeure partie des points de ce sous-ensemble à condition de ne pas être trop loin de la frontière. C'est ce qu'on appelle les SVM linéaire à marge légère (soft margin). De manière heuristique, comme en pratique, bien souvent, on ne peut pas avoir y_i(wx_i +b)-1  ≥ 0 pour tout i∈{ 1, ⋯, n }, on relâche en introduisant des variables positives ξ telle que 

    {[ ω^ Tx_i+b ≥ +1-ξ_i   lorsque   y_i=+1; ω^ Tx_i+b ≤ -1+ξ_i   lorsque   y_i=-1;                                       ].

avec ξ_i≥ 0. On a une erreur de classification si ξ_i>1, et on va alors introduire une pénalité, un coût à payer pour chaque erreur commise. On cherche alors à résoudre un problème quadratique

    min{1/2ω^ Tω+C1^T1_ξ>1}

sous la contrainte (<ref>), qui pourra être résolu de manière numérique très efficacement par descente de coordonnées (décrit auparavant).

S'il n'est pas possible de séparer les points, une autre astuce possible consiste à les transformer dans une dimension supérieure, de sorte que les données deviennent alors linéairement séparables. Trouver la bonne transformation qui sépare les données est toutefois très difficile. Cependant, il existe une astuce mathématique pour résoudre ce problème avec élégance, en définisant les transformations T (·) et les produits scalaires via un noyau
K (x _1, x _2) = ⟨ T (x _1),T (x _2)⟩.
L'un des choix les plus courants pour une fonction de noyau est la fonction de base radiale (noyau gaussien)
K (x _1, x _2) = exp (-‖x _1- x _2) ‖ ^ 2 ).
Il n'existe néanmoins pas de règles à ce jour permettant de choisir le « meilleur » noyau. Comme mentionné au début de la section précédante, cette technique est basé sur de la minimisation de distance, et il n'a aucune prévision de la probabilité d'être positif ou négatif (mais une interprétation probabiliste est néanmoins possible, comme le montre Grandvalet, par exemple).





 §.§ Arbres, Bagging et Forêts Aléatoires


Les arbres de classification ont été introduits dans Quinlan mais c'est surtout Breiman qui a assuré la popularité de l'algorithme. On parle de modèle CART pour « Classification And Regression Tree ». 
L'idée est de diviser consécutivement (par une notion de branchement) les données d'entrée jusqu'à ce qu'un critère d'affectation (par rapport à la variable cible) soit atteint, selon une règle prédéfinie. 

L'intuition de la construction des arbres de classification est la suivante. L'entropie H (x) est associée à la quantité de désordre dans les données x par rapport aux modalités prises par la variable de classification y, et chaque partition vise à réduire ce désordre. L'interprétation probabiliste est de créer les groupes les plus homogènes possible, en réduisant la variance par groupe (variance intra), ou de manière équivalente en créant deux groupes aussi différents que possible, en augmentant la variance entre les groupe (variance inter).
À chaque étape, nous choisissons la partition qui donne la plus forte réduction de désordre (ou de variance). L'arbre de décision complet se développe en répétant cette procédure sur tous les sous-groupes, où chaque étape k aboutit à une nouvelle partition en 2 branches, qui subdivise notre ensemble de données en 2. Enfin, on décide quand mettre fin à cette constitution de nouvelles branches, en procédant à des affectations finales (nœuds dits foliaires). Il existe plusieurs options pour mettre fin à cette croissance. L'une est de construire un arbre jusqu'à ce que toutes les feuilles soient pures, c'est à dire composées d'une seule observation. Une autre option est de définir une règle d'arrêt liée à la taille, ou à la décomposition, des feuilles. Les exemples de règles d'arrêt peuvent être d'une taille minimale (au moins 5 éléments par feuille), ou une entropie minimale. On parlera alors d'élagage de l'arbre: on laisse l'arbre grossir, puis on coupe certaines branches a posteriori (ce qui est différent de l'introduction d'un critère d'arrêt a priori au processus de croissance de l'arbre - par exemple en imposant une taille minimale aux feuilles, ou d'autres critères discutés dans Breiman).

À un nœud donné, constitué de n_0 observations (x_i,y_i) avec i∈ℐ_0, on va couper en deux branches (une à gauche et une à droite), partitionnant ainsi ℐ_0 en ℐ_g et ℐ_d. Soit I le critère d'intérêt, comme l'entropie du nœud (ou plutôt du nœud vu en tant que feuille):

    I(y_0)=-n_0 p_0log p_0    où  p_0=1/n_0∑_i∈ℐ_0y_i,

ou la variance du nœud:

    I(y_0)=n_0 p_0(1- p_0)    où  p_0=1/n_0∑_i∈ℐ_0y_i,

ce dernier étant également l'indice d'impureté de Gini. 



On partitionnera entre la branche gauche et la branche droite si le gain I(y_0)-[I(y_g)+I(y_d)] 
est suffisamment important. Lors de la construction des arbres, on va chercher la partition qui donne le gain le plus important possible. Ce problème combinatoire étant complexe, le critère suggéré par Breiman est de considérer un découpage suivant une des variables, avec ℐ_g={ i∈ℐ_0:x_k,i<s} et ℐ_d={ i∈ℐ_0:x_k,i>s}, pour une variable k et un seuil s (si la variable est continue, sinon on considère des regroupements de modalités pour des variables qualitatives). 

Les arbres de décision ainsi décrits sont simples à obtenir et faciles à interpréter (comme le montre la Figure <ref> sur les données du Titanic[Ce jeu de données, contenant des informations sur tous les passagers (et membres d'équipage) du Titanic, dont la variable y indiquant si la personne a survécu a été abondamment utilisé pour illustrer les techniques de classification, voir https://www.kaggle.com/c/titanic/data.]), mais ils sont peu robustes, et leur pouvoir prédictif est souvent très faible, en particulier si l'arbre est très profond.
Une idée naturelle est de développer un ensemble de modèles d'arbres à peu près indépendants, qui prédisent conjointement mieux qu'un modèle d'arbre unique. On va utiliser le bootstrap, en tirant (avec remise) n observations parmi {(_i,y_i)}. Ces ensembles d'arbres - naturellement appelés « forêts » - une fois agrégés donnent souvent de bien meilleurs résultats que les arbres isolés, mais elles sont difficiles à interpréter. Ces techniques ressemblent toutefois beaucoup à ce qui est fait lorsque l'on utilise les techniques de bootstrap en régression (par exemple pour construire des tubes de confiance dans une régression fonctionnelle).

Le principe du « bagging », pour « bootstrap aggregating »,  consiste à générer des échantillons aléatoires, en tirant avec remise dans l'échantillon d'origine, comme pour le bootstrap. Chaque échantillon ainsi généré  permet d'estimer un nouvel arbre de classification, formant ainsi une forêt d'arbres. C'est l'aggrégation de tous ces arbres qui conduit à la prévision. Le résultat global est moins sensible à l'échantillon initial et donne souvent de meilleurs résultats de prévision.

Les forêts aléatoires, ou « random forests » reposent sur le même principe que le « bagging », mais en plus, lors de la construction d'un arbre de classification, à chaque branche, un  sous-ensemble de m covariables est tiré aléatoirement. Autrement dit, chaque branche d'un arbre ne s'appuie pas sur le même ensemble de covariables. Cela permet d'amplifier la variabilité entre les différents arbres et d'obtenir, au final, une forêt composée d'arbres moins corrélés les uns aux autres. 



 §.§ Sélection de modèle de classification


Étant donné un modèle m(·) approchant [Y|=], et un seuil s∈[0,1], posons

    y^(s)=1[m()>s]
    =
    {[  1  si m()>s; 0  si m()≤ s ].

La matrice de confusion est alors le tableau de contingence associé aux comptages N=[N_u,v] avec

    N_u,v^(s)=∑_i=1^n 1(y^(s)_i=u,y_j=v)

pour (u,v)∈{0,1}. La Table <ref> présente un tel tableau, avec le nom de chacun des éléments : TP (true positive) sont les vrais positifs, correspondant aux 1 prédit en 1, 
TN (true negative) sont les vrais négatifs, correspondant aux 0 prédit en 0, FP (false positive) sont les faux positifs, correspondant aux 0 prédit en 1, et enfin FN (false negative) sont les faux négatifs, correspondant aux 1 prédit en 0).



Plusieurs quantités sont dérivées de ce tableau. La sensibilité correspond à la probabilité de prédire 1 dans la population des 1, ou taux de vrais positifs. La spécificité est  la probabilité de prédire 0 dans la population des 0 ou taux de vrais négatifs. On s'intéressera toutefois davantage au taux de faux négatifs,  c'est à dire la probabilité de prédire 1 dans la population des 0. La représentation de ces deux valeurs lorsque s varie donne la courbe ROC (« receiver operating characteristic ») :

    ROC_s=(FP_s/FP_s+TN_s,TP_s/TP_s+FN_s)=( sensibility_s,1- specificity_s)  pour s∈[0,1].

Une telle courbe est présentée dans la partie suivante, sur des données réelles.

Les deux grandeurs intensivement utilisées en machine learning sont l'indice κ, qui compare la précision observée avec celle espérée, avec un modèle aléatoire (tel que décrit dans LandisKoch) et l'AUC correspondant à l'aire sous la courbe ROC. Pour le premier indice, une fois choisi s, notons N^⊥ le tableau de contingence correspond aux cas indépendants (défini à partir de N dans le test d'indépendance du chi-deux).
On pose alors

    précision totale=TP+TN/n

alors que

    précision aléatoire=[TN+FP]·[TP+FN] 
    +[TP+FP]·[TN+FN]/n^2

On peut alors définir

    κ=précision totale-précision aléatoire/1-précision aléatoire

Classiquement s sera fixé égal à 0.5, comme dans une classification bayésienne naïve, mais d'autres valeurs peuvent être retenues, en particulier si les deux erreurs ne sont pas symmétriques (nous reviendrons sur ce point dans un exemple par la suite).

Il existe des compromis entre des modèles simples et complexes mesurés par leur nombre de paramètres (ou plus généralement les degrés de liberté) en matière de performance et de coût. Les modèles simples sont généralement plus faciles à calculer, mais peuvent conduire à des ajustements plus mauvais (avec un biais élevé par exemple). Au contraire, les modèles complexes peuvent fournir des ajustements plus précis, mais risquent d'être coûteux en termes de calcul. En outre, ils peuvent surpasser les données ou avoir une grande variance et, tout autant que des modèles trop simples, ont de grandes erreurs de test. Comme nous l'avons rappelé auparavant, dans l'apprentissage machine, la complexité optimale du modèle est déterminée en utilisant le compromis de biais-variance.



 §.§ De la classification à la régression


Comme nous l'avons rappelé en introduction, historiquement, les méthodes de machine learning se sont orientées autour des problèmes de classification (avec éventuellement plus de 2 modalités[Par exemple dans le cas de  reconnaissance de lettres ou de chiffres]), et assez peu dans le cas où la variable d'intérêt y est continue. Néanmoins, il est possible d'adapter quelques techniques, comme les arbres et les forêts aléatoires, le boosting, ou les réseaux de neurones.

Pour les arbres de régression, Morganetal ont proposé la méthode AID, basée sur la formule de décomposition de la variance de l'équation (<ref>), avec un algorithme proche de celui de la méthode CART décrite auparavant. Dans le contexte de la classification, on calculait, à chaque nœud (dans le cas de l'indice d'impureté de Gini) en sommant sur la feuille de gauche {x_k,i<s} et celle de droite {x_k,i>s}

    I=∑_i:x_k,i<sy_g(1-y_g)
    +∑_i:x_k,i>sy_d(1-y_d)

où y_g et y_d désignent les fréquences de 1 dans la feuille de gauche et de droite, respectivement. Dans le cas d'un arbre de régression, on utilisera

    I=∑_i:x_k,i<s(y_i-y_g)^2
    +∑_i:x_k,i>s(y_i-y_d)^2

qui va correspondre à la somme (pondérée) des variances intra. Le partage optimal sera celui qui aura le plus de variance intra (on veut les feuilles les plus homogènes possibles) ou de manière équivalente, on veut maximiser la variance intra.

Dans le contexte des forêts aléatoires, on utilise souvent un critère majoritaire en classification (la classe prédite sera la classe majoritaire dans une feuille), alors que pour la régression, on utilise la moyenne des prédictions, sur tous les arbres.




Dans la partie précédente, nous avons présenté la dimension « apprentissage » du machine learning en présentant le boosting. Dans un contexte de régression (variable y continue), l'idée est de créer une succession de modèles en écrivant l'équation (<ref>) sous la forme :

    m^(k)(x)=m^(k-1)(x)+α_k h∈ℋargmin{∑_i=1^n (y_i,m^(k-1)(x)+h(x))^2}

où α_k est un paramètre de « shrinkage », où le second terme correspond à un arbre de régression, sur les résidus, y_i-m^(k-1)(x_i).

Mais il existe d'autres techniques permettant d'apprendre de manière séquentielle. Dans un modèle additif (gam) on va chercher une écriture de la forme

    m(x)=∑_j=1^p m_j(x_j)=m_1(x_1)+⋯+m_p(x_p)

L'idée de la poursuite de projection repose sur une décomposition non pas sur les variables explicatives, mais sur des combinaisons linéaires. On va ainsi considérer un modèle

    m(x)=∑_j=1^k g_j(ω_jx)=g_1(ω_1x)+⋯+g_k(ω_kx).

Tout comme les modèles additifs, les fonctions g_1,⋯,g_k sont à estimer, tout comme les directions ω_1,⋯,ω_k. Cette écriture est relativement générale, et permet de tenir compte d'intéractions et d'effets croisés (ce que nous ne pouvions pas faire avec les modèles additifs qui ne tiennent compte que de non-linéarités). Par exemple en dimension 2, un effet multiplicatif m(x_1,x_2)=x_1· x_2 s'écrit

    m(x_1,x_2)=x_1· x_2=(x_1+x_2)^2/4-(x_1-x_2)^2/4

autrement dit g_1(x)=x^2/4, g_1(x)=-x^2/4, ω_1=(1,1) et ω_1=(1,-1). Dans la version simple, avec k=1, avec une fonction de perte quadratique, on peut utiliser un développement de Taylor pour approcher [y_i-g(ωx_i)]^2, et construire classiquement un algorithme itératif. Si on dispose d'une valeur initiale ω_0, notons que

    ∑_i=1^n[y_i-g(ωx_i)]^2
    ≈∑_i=1^n g'(ω_0x_i)^2[ωx_i+y_i-g(ω_0x_i)/g'(ω_0x_i)-
    ωx_i]^2

qui correspondrait à l'approximation dans les modèles linéaires généralisés sur la fonction g(·) qui était la fonction de lien (supposée connue). On reconnait  un problème de moindres carrés pondérés. La difficulté ici est que les fonctions g_j(·) sont inconnues.



§ APPLICATIONS


Les données massives ont rendu nécessaire le développement de techniques d'estimation permettant de pallier les limites des modèles paramétriques, jugés trop restrictifs, et  des modèles non-paramétriques classiques, dont l'estimation peut être difficile en présence d'un nombre élevé de variables. L'apprentissage statistique, ou apprentissage machine, propose de nouvelles méthodes d'estimation non-paramétriques, performantes dans un cadre général et en présence d'un grand nombre de variables.[Entre autres, voir HastieEtal et James.] Toutefois, l'obtention d'une plus grande flexibilité s'obtient au prix d'un manque d'interprétation qui peut être important.

En pratique, une question importante est de savoir quel est le meilleur modèle ? La réponse à cette question dépend du problème sous-jacent. Si la relation entre les variables est correctement approximée par un modèle linéaire, un modèle paramétrique correctement spécifié devrait être performant. Par contre, si le modèle paramétrique n'est pas correctement spécifié, car la relation est fortement non-linéaire et/ou fait intervenir des effets croisés non-négligeables, alors les méthodes statistiques issues du machine learning devraient être plus performantes.

La bonne spécification d'un modèle de régression est une hypothèse souvent posée, elle est rarement vérifiée et justifiée. Dans les applications qui suivent, nous montrons comment les méthodes statistiques issues du machine learning peuvent être utilisées pour justifier la bonne spécification d'un modèle de régression paramétrique, ou pour détecter une mauvaise spécification. Des applications en classification sont présentées dans un premier temps, sections <ref>, <ref> et <ref>. D'autres applications sont ensuite présentées dans le contexte de régression classique, sections <ref> et <ref>.



 §.§ Les ventes de sièges auto pour enfants (classification)



Nous reprenons ici un exemple utilisé dans James. Le jeu de données contient les ventes de sièges auto pour enfants dans 400 magasins (Sales), ainsi que plusieurs variables, dont la qualité de présentation en rayonnage (Shelveloc, égal à « mauvais », « moyen », « bon ») et le prix (Price).[C'est le jeu de données Carseats de la bibliothèque ISLR.] Une variable dépendante binaire est artificiellement crée, pour qualifier une forte vente ou non (High=« oui » si >8 et à « non » sinon). Dans cette application, on cherche à évaluer les déterminants d'un bon niveau de vente.

Dans un premier temps, on considère un modèle de régression linéaire latent:

    y^⋆ = γ +  + ε,     ε∼ G(0,1),

où  est composé de k  variables explicatives,  est un vecteur de k paramètres inconnus et ε est un terme d'erreur i.i.d. avec une fonction de répartition G d'espérance nulle et de variance égale à un. La variable dépendante y^⋆ n'est pas observé, mais seulement y, avec:

    y=
    
    1    si  y^⋆ > ξ, 
    
    0    si  y^⋆≤ξ.

On peut alors exprimer la probabilité d'avoir y égal à 1, comme suit :

    ℙ(Y=1)=G(β_0+)

où β_0=γ-ξ.[ℙ[Y=1]=ℙ[Y^⋆>ξ]=ℙ[γ +  + ε>ξ]=ℙ[ ε>ξ-γ - ]=ℙ[ε<γ-ξ + ]. En posant γ-ξ=β_0, on obtient ℙ[Y=1]=G(β_0+). En général, on suppose que le terme d'erreur est de variance σ^2, auquel cas les paramètres du modèle (<ref>) deviennent β_0/σ et /σ, ce qui veut dire que les paramètres du modèle latent (<ref>) ne sont pas identifiables, ils sont estimés à un paramètre d'échelle près.] L'estimation de ce modèle se fait par maximum de vraisemblance en sélectionnant a priori une loi paramétrique G. Si on suppose que G est la loi Normale, c'est un modèle probit, si on suppose que G est la loi logistique, c'est un modèle logit. Dans un  modèle  logit/probit, il y a deux sources potentielles de mauvaise spécification :


  (i) la relation linéaire β_0+ est mal spécifiée

  (ii) la loi paramétrique utilisée G n'est pas la bonne

En cas de mauvaise spécification, de l'une ou l'autre sorte, l'estimation n'est plus valide. Le modèle le plus flexible est le suivant :

    ℙ[Y=1|X=x]=G(h(x))

où h est une fonction inconnue et G une fonction de répartition inconnue. Les modèles de bagging, de forêt aléatoire et de boosting permettent d'estimer ce modèle  général sans faire de choix à priori sur la fonction h  et sur la distribution G. L'estimation du modèle logit/probit est néanmoins plus performante si h et G sont correctement spécifiés.






Nous estimons le modèle (<ref>) avec la loi logistique pour G, et le modèle (<ref>) avec les méthodes de bagging, de forêt aléatoire et de boosting. Nous faisons une analyse de validation croisée par 10 blocs. Les probabilités individuelles des données out-of-sample, c'est à-dire de chacun des blocs non-utilisée pour l'estimation, sont utilisées pour évaluer la qualité de la classification. 

La figure <ref> présente la courbe ROC, ainsi que l'aire sous la courbe (AUC), pour les estimations logit, bagging, random forest et boosting. La courbe ROC est un graphique qui représente simultanément la qualité de la prévision dans les deux classes, pour des valeurs différentes du seuil utilisé pour classer les individus (on parle de « cutoff »). Une manière naturelle de classer les individus consiste à les attribuer dans la classe pour laquelle ils ont la plus grande probabilité estimée. Dans le cas d'une variable binaire, cela revient à prédire la classe d'appartenance pour laquelle la probabilité estimée est supérieure à 0.5. Mais un autre seuil pourrait être utilisé. Par exemple, dans la figure <ref>, un point de la courbe ROC du modèle logit indique qu'en  prenant un seuil égal à 0.5, la réponse "non" est correctement prédite à 90.7% (specificity), et la réponse « oui »  à  86% (sensitivity). Un autre  point indique qu'en prenant un seuil égal à 0.285, la réponse « non » est correctement prédite à 86% (specificity), et la réponse « oui »  à  92.7% (sensitivity). Comme décrit auparavant, un modèle de classification idéal aurait une courbe ROC de la forme Γ. Autrement dit, le meilleur modèle est celui dont la courbe est au-dessus des autres. Un critère souvent utilisé pour sélectionner le meilleur modèle est celui dont l'aire sous la courbe ROC est la plus grande (AUC). L'avantage d'un tel critère est qu'il est simple à comparer et qu'il ne dépend pas du choix du seuil de classification.

Dans notre exemple, la courbe ROC du modèle logit domine les autres courbes, et son aire sous la courbe est la plus grande (AUC=0.9544). Ces résultats indiquent que ce modèle  fournit les meilleures prévisions de classification. N'étant dominé par aucun autre modèle, ce constat suggère que le modèle linéaire logit est correctement spécifié et qu'il n'est pas utile d'utiliser un modèle plus général et plus complexe.




 §.§ L'achat d'une assurance caravane (classification)



Nous reprenons à nouveau un exemple utilisé dans James. Le jeu de données contient 85 variables sur les charactéristiques démographiques de 5822 individus.[C'est le jeu de données Caravan de la bibliothèque ISLR sous R.]  La variable dépendante (Purchase) indique si l'individu a acheté une assurance caravane, c'est une variable binaire, égale à « oui » ou « non ». Dans le jeu de données, seulement 6% des individus ont pris une telle assurance. Les classes sont donc fortement désequilibrées.






Nous estimons le modèle (<ref>) avec la loi logistique et le modèle (<ref>) avec les méthodes bagging, forêt aléatoire et boosting (les paramètres de « tuning » sont ceux de James, n.trees=1000 et shrinkage=0.01).  Nous faisons une analyse de validation croisée par 10 blocs. Les probabilités individuelles des données out-of-sample, c'est à-dire de chacun des blocs non-utilisée pour l'estimation, sont utilisées pour évaluer la qualité de la classification. 


La figure <ref> présente la courbe ROC, ainsi que l'aire sous la courbe (AUC), pour les estimations logit, bagging, random forest et boosting. La courbe du modèle boosting domine les autres courbes, son aire sous la courbe est la plus grande (AUC=0.7691). Ces résultats indiquent que le boosting  fournit les meilleures prévisions de classification. Notons que, comparées à l'exemple précédent,  les courbes sont assez éloignées de la forme en coude,  ce qui suggère que la classification ne sera pas aussi bonne.

Il faut faire attention aux résultats d'une classification standard, c'est-à-dire avec un seuil de classification égal à 0.5, qui est souvent pris par défaut dans les logiciels (la prédiction de la réponse de l'individu i est « non » 1 si la probabilité estimée qu'il réponde « non » est supérieure à 0.5, sinon c'est « oui »). La partie gauche du tableau <ref> présente les taux de classifications correctes avec ce seuil (0.5 cutoff), pour les différentes méthodes. Avec le meilleur modèle et le seuil standard (boosting et seuil à 0.5), les réponses « non » sont correctes à 99.87% (spécificité, specificity) et les réponses « oui » sont toutes fausses (sensitivité, sensitivity). Autrement dit, cela équivaut à utiliser un modèle qui prédit que personne n'achète d'assurance caravane. Sélectionner un tel modèle est absurde pour l'analyste, qui est surtout intéressé par les 6% des individus qui en ont pris une. Ce résultat s'explique par la présence de classes fortement déséquilibrées. En effet, dans notre exemple, en prévoyant que personne n'achète d'assurance, on  fait « seulement » 6% d'erreur. Mais ce sont des erreurs qui conduisent à ne rien expliquer.


	

Plusieurs méthodes peuvent être utiles pour pallier à ce problème, lié aux classes fortement désequilibrées (pour plus d'informations, voir KuhnJohnson, chapitre 16). Une solution simple consiste à utiliser un seuil de classification différent. La courbe ROC présente les résultats en fonction de plusieurs seuils de classification, où la classification parfaite est illustrée par le couple (specificity, sensitivity)=(1,1), c'est à-dire par le coin supérieur gauche dans le graphique. Aussi, on choisit comme seuil de classification optimal celui qui correspond au point de la courbe ROC qui est le plus proche du point (1,1), ou du coin supérieur gauche. La partie droite du tableau <ref> présente les taux de classifications correctes avec les seuils optimaux (optimal cutoff), pour les différentes méthodes (les seuils optimaux des méthodes logit, bagging, forêt aléatoire et boosting sont, respectivement, égaux à 0.0655 , 0.0365 , 0.0395, 0.0596). Avec le boosting et un seuil optimal, les réponses "non" sont correctes à 68.6% (specificity) et les réponses «  oui » à 73.85%  (sensitivity). L'objet de l'analyse étant de prévoir correctement les individus susceptibles d'acheter une assurance caravane (classe "oui"), et les distinguer suffisamment des autres (classe «  non »), le choix du seuil optimal est beaucoup plus performant que le seuil standard 0.5. Notons qu'avec un modèle logit et un seuil optimal, le taux de classifications correctes de la classe "non" est de 72.78%, celui de la classe "oui" est de 63.51%. Par rapport au boosting, le logit prédit un peu mieux la classe "non", mais nettement moins bien la classe «  oui ».



 §.§ Les défauts de remboursement de crédits particuliers (classification)




Considérons la base allemande de crédits particuliers, utilisée dans Nisbet et Tuffery, avec 1000 observations, et 19 variables explicatives, dont 12 qualitatives c'est à dire, en les disjonctant (en créant une variable indicatrice pour chaque modalité), 48 variables explicatives potentielles.

Une question récurrente en modélisation est de savoir quelles sont les variables qui mériteraient d'être utilisées. La réponse la plus naturelle pour un économètre pourrait être une méthode de type stepwise (parcourir toutes les combinaisons possibles de variables étant a priori un problème trop complexe en grande dimension). La suite des variables dans une approche forward est présentée dans la première colonne du tableau <ref>. Une approche mentionnée avant qui peut être utile est le Lasso, en pénalisant convenablement la norme ℓ_1 du vecteur de paramètres . On peut ainsi, séquentiellement, trouver les valeurs du paramètre de pénalisation λ, qui permet d'avoir une variable explicative supplémentaire, non nulle. Ces variables sont présentées dans la dernière colonne. On note que les deux premières variables considérées comme non nulles (pour un λ assez grand) sont les deux premières à ressortir lors d'une procédure forward. Enfin, une dernière méthode a été proposée par  Breiman2001, en utilisant tous les arbres créé lors de la construction d'une forêt aléatoire : l'importance de la variable x_k dans une forêt de T arbres est donnée par:

    Importance(x_k)=1/T∑_t=1^n
    ∑_j∈ N_t,k p_t(j)Δℐ(j)

où N_t,k désigne l'ensemble des nœuds de l'arbre t utilisant la variable x_k comme variable de séparation, p_t(j) désigne la proportion des observations au nœud j, et Δ(j) est la variation d'indice au nœud j (entre le nœud précédant, la feuille de gauche et celle de droite). Dans la colonne centrale du tableau <ref> sont présentées les variables par ordre d'importance décroissante, lorsque l'indice utilisé est l'indice d'impureté de Gini.



Avec l'approche stepwise et l'approche lasso, on reste sur des modèles logistiques linéaires. Dans le cas des forêtes aléatoires (et des arbres), des intéractions entre variables peuvent être prises en compte, lorsque 2 variables sont présentes. Par exemple la variable residence_since est présente très haut parmi les variables prédictives (troisième variable la plus importante).












































 §.§ Les déterminants des salaires (régression)
 


Afin d'expliquer les salaires (individuels) en fonction du niveau d'étude, de l'expérience de la personne, et son genre, il est classique d'utiliser l'équation de salaire de Mincer - décrite dans Mincer - tel que le rappelle Lemieux:

    log () = β_0+β_1  +β_2  +β_3  ^2+β_4  +ε

où ed est le niveau d'études, ex l'expérience professionnelle et fe une variable indicatrice, égale à 1 si l'individu est une femme et à 0 sinon. D'après la théorie du capital humain, le salaire espéré augmente avec l'expérience, de moins en moins vite, pour atteindre un maximum avant de diminuer. L'introduction du carré de exp permet de prendre en compte une telle relation. La présence de la variable  fe permet quand à elle de mesurer une éventuelle discrimination salariale entre les hommes et les femmes.

Le modèle (<ref>) impose une relation linéaire entre le salaire  et le niveau d'étude, et une relation quadratique entre le salaire et l'expérience professionnelle. Ces relations peuvent paraître trop restrictives. Plusieurs études montrent notamment que le salaire ne diminue pas après un certain age, et qu'une relation quadratique ou un polynôme de degré plus élevé est plus adapté (comme décrit dans MurphyWelch et  Bazen).

Le modèle (<ref>) impose également que la différence salariale entre les hommes et les femmes est indépendente du niveau d'étude et de l'expérience. Il est trop restrictif si, par exemple, on suspecte que l'écart de salaire moyen entre les hommes et les femmes est faible pour les postes non-qualifiés et fort pour les postes qualifiés, ou faible en début de carrière et fort en fin de carrière (effets d'intéractions). 

Le modèle le plus flexible est le modèle entièrement non-paramétrique :

    log () = m() +ε

où m(·) est une fonction quelconque. Il a l'avantage de pouvoir tenir compte de relations non-linéaires quelconques et d'intéractions complexes entre les variables. Mais, sa grande flexibilité se fait au détriment d'une interprétation plus difficile du modèle. En effet, il faudrait un graphique en 4-dimensions pour représenter la fonction m. Une solution consiste à représenter la fonction m en 3 dimensions, en fixant la valeur de l'une des variables, mais la fonction représentée peut  être très différente avec une valeur fixée différente.







Nous utilisons les données d'une enquête de l'US Census Bureau daté de mai 1985, issues de l'ouvrage de Berndt et disponibles sour R.[C'est le jeu de données CPS1985 de la bibliothèque AER.] Nous estimons les deux modèles et utilisons une analyse de validation croisées par 10 blocs pour sélectionner la meilleure approche. Le modèle paramétrique (<ref>) est estimé par Moindres Carrés Ordinaires (ols). Le modèle entièrement non-paramétrique (<ref>) est estimé par la méthode des splines,  car il en comprend peu de variables, ainsi que par les méthodes bagging, random forest et boosting.

Le tableau <ref> présente les résultats de la  validation croisée en 10 blocs (10-fold cross-validation). Le meilleur modèle est celui qui minimise le critère 
ℛ^10- CV. Les résultats montrent que le modèle (<ref>) est au moins aussi performant que le modèle (<ref>), ce qui suggère que le modèle paramétrique (<ref>) est correctement spécifié.  





 §.§ Les déterminants des prix des logements à Boston (régression)




Nous reprenons ici l'un des exemples utilisé dans James, dont les données sont disponibles sous R. Le jeu de données contient les valeurs médianes des prix des maisons (medv) dans n=506 quartiers autour de Boston, ainsi que 13 autres variables, dont le nombre moyen de pièces par maison (rm), l'age moyen des maisons (age) et le pourcentage de ménages dont la catégorie socio-professionnelle est peu élevée (lstat).[C'est le jeu de données  de  la librairie . Pour une description complète des données, voir: https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/Boston.html.] 

Considérons le modèle de régression linéaire suivant :

    = α +  +ε

où =[] est un vecteur en dimension 13 et  est un vecteur de 13 paramètres. Ce modèle spécifie une relation linéaire entre la valeur des maisons et chacune des variable explicatives. 


Le modèle le plus flexible est le modèle entièrement non-paramétrique :

    = m( ) +ε.

L'estimation de ce modèle avec les méthodes du noyau ou les splines peut être problématique, car le nombre de variables est relativement élevé (il y a ici 13 variables), ou au moins trop élevé pour envisager estimer une surface en dimension 13. 
Nous estimons les deux modèles et utilisons une analyse de validation croisée par 10-blocs pour sélectionner la meilleure approche. Le modèle paramétrique (<ref>) est estimé par Moindres Carrés Ordinaires (ols) et le modèle entièrement non-paramétrique (<ref>) est estimé par trois méthodes différentes: bagging, forêt aléatoire et boosting (nous utilisons ici les valeurs par défaut utilisées dans James, pp. 328-331).




Le tableau <ref> présente les résultats de la  validation croisée en 10 blocs (10-fold cross-validation). La première ligne (in-sample) présente la qualité de l'ajustement des modèles en utilisant seulement les données d'apprentissage, c'est-à-dire celles qui ont servi à estimer le modèle, pour calculer le mse. La deuxième ligne (out-of-sample) présente la qualité de l'ajustement en utilisant d'autres données que celles ayant servies à estimer le modèle, pour calculer l'erreur quadratique.  À partir  des résultats in-sample, les méthodes de bagging et de random forest paraissent incroyablement plus performantes que l'estimation ols du modèle linéaire (<ref>), le critère 
ℛ^10- CV passant de 21.782 à 1.867 et 1.849.
Les résultats out-of-sample vont dans le même sens, mais la différence est  moins importante, le critère 
ℛ^10- CV passant de 24.082 à 9.59 et 9.407. Ces résultats illustrent un phénomène classique des méthodes non-linéaires, comme le bagging et la forêt aléatoire, qui peuvent être très performantes pour prédire les données utilisées pour l'estimation, mais  moins performantes pour prédire des données hors-échantillon. C'est pourquoi la sélection de la meilleure estimation est habituellement basée sur une analyse out-of-sample, telle que présentée dans la deuxième ligne.

La différence entre l'estimation du modèle linéaire (<ref>) et du modèle entièrement non-paramétrique (<ref>) est importante (24.082 vs 9.590, 9.407 et 11.789). Un tel écart suggère que le modèle linéaire est mal spécifié, et que des relations non-linéaire et/ou des effets d'intéractions sont présentes dans la relation entre le prix des logements, medv, et les variables explicatives . Ce résultat nous conduit à  chercher une meilleure spécification paramétrique.


À partir du modèle paramétrique  (<ref>), et afin de prendre en compte d'éventuelles non-linéarités, le modèle additif généralisé (gam) suivant peut être considéré :

    = m_1 (x_1)+m_2 (x_2)+…+m_13 (x_13) +ε,

où m_1,  m_2, … m_13 sont des fonctions inconnues. L'avantage de ce modèle est qu'il permet de considérer n'importe quelle relation non-linéaire entre la variable dépendante et chacune des variables explicatives. De plus, il ne souffre pas du problème du fléau de la dimension, car chacune des fonction est de dimension 1, et il est facilement interprétable. Toutefois, il ne prend pas en compte d'éventuels effets d'intéractions. 

L'estimation du modèle additif généralisé (<ref>) par la méthode des splines, dans le cadre d'une analyse de validation croisée par 10-blocs, donne une valeur 
ℛ^10- CV=13.643. Par rapport au modèle paramétrique (<ref>), il y a un gain important (13.643 vs. 24.082). Mais la différence avec le modèle entièrement non-paramétrique (<ref>) reste conséquente (13.643 vs 9.590, 9.407, 11.789). Une telle différence suggère que la prise en compte de relations individuelles pouvant être fortement non-linéaires n'est pas suffisante, et que des effets d'intéractions entre les variables sont présents. Nous pourrions inclure dans le modèle les variables d'intéractions les plus simples entre toutes les paires de variables (x_i× x_j), mais cela impliquerait de rajouter un très grand nombre de variables au modèle initial (78 dans notre cas), qui ne serait pas sans conséquence sur la qualité de l'estimation du modèle. Quoi qu'il en soit, nous pouvons dire pour le moment que le modèle linéaire est mal spécifié et qu'il existe des effets d'intéractions pouvant être forts dans la relation entre medv et X, l'identification de tels effets restant délicat. 



Afin d'aller plus loin, les outils développés en apprentissage statistique peuvent être à nouveau d'un grand recours. Par exemple, l'estimation random forest s'accompagne de mesures de l'importance de chacune des variables dans l'estimation du modèle (décrit dans la section précédante). Le tableau <ref> présente ces mesures dans le cadre du modèle (<ref>), estimé sur  l'échantillon complet. Les résultats suggèrent que les variables rm et lstat sont les variables les plus importantes pour expliquer les variations des prix des logements medv.
Ce constat nous conduit à enrichir la relation initiale, en rajoutant les intéractions liées à ces deux variables seulement, qui sont les plus importantes. 




Nous estimons le modèle additif généralisé incluant les variables d'intéractions, sur l'échantillon complet:

    = m_1(x_1)+ ⋯ + m_13(x_13)+ (:x)  γ + (:x)  δ +ε,

où (:x) représente les variables d'intéractions de rm avec toutes les autres variables de x et (:x) représente les variables d'intéractions de lstat avec toutes les autres variables de x.[On a 
(:x)=[rm×chas, rm×nox,  rm×age, rm×tax, rm×indus, 
rm×rad, rm×dis, rm×lstat, rm×crim, rm×black, rm×zn, rm×ptratio] et 
(:x)=[lstat×chas, lstat×nox, lstat×age, lstat×tax, lstat×indus, 
lstat×rad, lstat×dis, lstat×crim, lstat×black, lstat×zn, lstat×ptratio].
]
L'analyse des résultats de cette estimation suggère que les fonctions m̂_i sont linéaires pour toutes les variables, sauf pour la variable dis, dont la relation estimée est présentée dans la figure <ref>.  Cette variable mesure la distance moyenne à cinq centres d'emplois de la région. L'effet semble diminuer plus rapidement avec la distance, lorsque celle-ci n'est pas très élevée. Au delà d'une certaine distance (au delà de 2, en log), l'effet est réduit, il continue à diminuer mais plus doucement. Cette relation non-linéaire peut  être approchée par une régression linéaire par morceaux en considérant un nœud.

Finalement, l'analyse précédente nous conduit à considérer  le modèle linéaire suivant:

    = α + + (-2)_+  θ + (: x)  γ + (:x) δ +ε

où (-2)_+ est égal à  la valeur de son argument si ce dernier est positif, et à 0 sinon. Par rapport au modèle linéaire initial, ce modèle inclut une relation linéaire par morceaux avec la variable dis, ainsi que des effets d'intéractions entre rm, lstat et chacune des autres variables de . 

Le tableau <ref> présente les résultats de la  validation croisée en 10 blocs (10-fold cross-validation) de l'estimation des modèles paramétriques (<ref>) et (<ref>),  estimés par Moindres Carrés Ordinaires (ols), et du modèle additif généralisé (<ref>) estimé par les splines. Il montre que l'ajout des variables d'intéractions et de la relation linéaire par morceaux dans le modèle (<ref>) donne des résultat beaucoup plus performants que le modèle initial (<ref>): le critère 
ℛ^10- CV est divisé par plus de deux, il passe de 24.082 à 11.759. En comparant ces résultats avec ceux du tableau  <ref>, on constate également que le modèle paramétrique (<ref>), estimé par ols, est aussi performant que le modèle général (<ref>) estimé par boosting (
ℛ^10- CV=11.789). La différence avec les méthodes bagging et forêt aléatoire n'est quant à elle pas très importante (
ℛ^10- CV=9.59, 9.407)

Finalement, les méthodes bagging, forêt aléatoire et boosting ont permis de mettre en évidence une mauvaise spécification du modèle paramétrique initial, puis de trouver un modèle paramétrique beaucoup plus performant, en prenant compte des effets de non-linéarités et d'intéractions appropriées.















































































































































































































































§ CONCLUSION


Si les « deux cultures » (ou les deux communautés) de l'économétrie et du machine learning se sont développées en parallèle, le nombre de passerelles entre les deux ne cesse d'augmenter. Alors que Varian présentait les apports importants de l'économétrie à la communauté du machine learning, nous avons tenté ici de présenter des concepts et des outils développés au fil du temps par ces derniers, qui pourraient être utiles aux économètres, dans un contexte d'explosion du volume de données. Si nous avons commencé par opposer ces deux mondes, c'est aussi pour mieux comprendre leurs forces et leurs faiblesses. Les fondements probabilistes de l'économétrie sont incontestablement sa force, avec non seulement une interprétabilité des modèles, mais aussi une quantification de l'incertitude. Néanmoins, nous l'avons vu à plusieurs reprises sur des données réelles, les performances prédictives des modèles de machine learning sont intéressantes, car elles permettent de mettre en avant une mauvaise spécification d'un modèle économétrique. De la même manière que les techniques non-paramétriques permettent d'avoir un point de référence pour juger de la pertinence d'un modèle paramétrique, les outils de machine learning permettent d'améliorer un modèle économétrique, en détectant un effet non-linéaire ou un effet croisé oublié.

Une illustration des interactions possibles entre les deux communautés se trouve par exemple dans BelloniBelloni2, dans un contexte de choix d'instrument dans une régression. Reprenant les données de  AngristKrueger sur un problème de réussite scolaire, ils montrent comment mettre en œuvre efficacement les techniques d'économétrie instrumentale quand on peut choisir parmi 1530 instruments disponibles (problème qui deviendra récurrent avec l'augmentation du volume de données).
Comme nous l'avons vu tout au long de cet article, même si les approches peuvent être fondamentalement différentes dans les deux communautés, bon nombre d'outils développés par la communauté du machine learning méritent d'être utilisés par les économètres.





99
[Ahamada & Flachaire (2011)]Ahamada Ahamada, I. & E. Flachaire (2011). Non-Parametric Econometrics. Oxford
University Press.
[Aigner et al. (1977)]Aigneretal Aigner, D., Lovell, C.A.J & Schmidt, P. (1977). Formulation and estimation of stochastic frontier production function models. Journal of Econometrics, 6, 21–37.
[Aldrich (2010)]Aldrich Aldrich, J. (2010). The Econometricians' Statisticians, 1895-1945. History of Political Economy, 42 111–154.
[Altman et al. (1994)]Altman Altman, E., Marco, G. & Varetto, F. (1994). Corporate distress diagnosis: Comparisons using linear discriminant analysis and neural networks (the Italian experience).Journal of Banking & Finance 18, 505–529.
[Angrist & Lavy (1999)]Angrist Angrist, J.D. & Lavy, V. (1999). Using Maimonides' Rule to Estimate the Effect of Class Size on Scholastic Achievement. Quarterly Journal of Economics, 114, 533–575. 
[Angrist, J.D. & Pischke, J.S. (2010)]AngristPischke10 Angrist, J.D. & Pischke, J.S. (2010). The Credibility Revolution in Empirical Economics: How Better Research Design Is Taking the Con out of Econometrics.  Journal of Economic Perspective, 24,  3–30.
[Angrist & Pischke (2015)]AngristPischke Angrist, J.D. & Pischke, J.S. (2015). Mastering Metrics. Princeton University Press.
[Angrist & Krueger (1991)]AngristKrueger Angrist, J.D. & Krueger, A.B. (1991). Does Compulsory School Attendance Affect Schooling and Earnings? Quarterly Journal of Economics, 106, 979–1014.
[Bottou (2010)]bottou Bottou, L. (2010) Large-Scale Machine Learning
with Stochastic Gradient Descent Proceedings of the 19th International Conference on Computational Statistics (COMPSTAT'2010), 177–187.
[Bajari et al. (2015)]Bajari Bajari, P., Nekipelov, D., Ryan, S.P. & Yang, M. 2015. Machine learning
methods for demand estimation. American Economic Review, 105 481–485.
[Bazen & Charni (2015)]Bazen Bazen, S. & K. Charni (2015). Do earnings really decline for older workers?
AMSE 2015-11 Discussion Paper, Aix-Marseille University.
[Bellman (1957)]Bellman Bellman, R.E. (1957). Dynamic programming. Princeton University Press.
[Belloni et al. (2010]Belloni Belloni, A., Chernozhukov, V. & Hansen, C. (2010). Inference Methods for High-Dimensional Sparse Econometric Models. Advances in Economics and Econometrics, 245–295
[, 2012)]Belloni2 Belloni, A., Chen, D., Chernozhukov, V. & Hansen, C. (2012). Sparse Models and Methods for Optimal Instruments With an Application to Eminent Domain. Econometrica, 80, 2369–2429. 
[Benjamini & Hochberg (1995)]Benjamini Benjamini, Y. & Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society, Series B, 57:289–300.
[Berger (1985)]Berger Berger, J.O. (1985). Statistical decision theory and Bayesian Analysis (2nd ed.). Springer-Verlag. 
[Berk (2008)]Berk Berk, R.A. (2008). Statistical Learning from a Regression Perspective. Springer Verlag.
[Berkson (1944)]Berkson1 Berkson, J. (1944). Applications of the logistic function to bioassay. Journal of the American
Statistical Association, 9, 357–365.
[Berkson (1951)]Berkson2 Berkson, J. (1951). Why I prefer logits to probits. Biometrics, 7 (4), 327–339.
[Bernardo & Smith (2000)]Bernardo Bernardo, J.M. & Smith, A.F.M. (2000). Bayesian Theory. John Wiley.
[Berndt (1990)]Berndt Berndt, E. R. (1990). The Practice of Econometrics: Classic and Contemporary. Addison Wesley.
[Bickel et al. (1997)]Bickel Bickel, P.J., Gotze, F. & van Zwet, W. (1997). Resampling fewer than n observations: gains, losses and remedies for losses. Statistica Sinica, 7, 1-31.
[Bishop (2006)]Bishop Bishop, C. (2006). Pattern Recognition and Machine Learning. Springer Verlag.
[Blanco et al. (2013)]Blanco Blanco, A. Pino-Mejias, M., Lara, J. & Rayo, S. (2013). Credit scoring models for the microfinance industry using neural networks: Evidence from peru. Expert Systems with Applications, 40, 356–364.
[Bliss (1934)]Bliss Bliss, C.I. (1934). The method of probits. Science, 79, 38–39.
[Breiman (2001a)]Breiman Breiman, L. (2001a). Statistical Modeling: The Two Cultures. Statistical Science, 16:3, 199–231.
[Bühlmann & van de Geer (2011)]Buhlmann Bühlmann, P. & van de Geer, S. (2011). Statistics for high-dimensional data: methods, theory and applications. Springer Verlag.
[Breiman (2001b)]Breiman2001 Breiman, L. (2001b). Random forests. Machine learning, 45:1, 5–32.

[Brown (1986)]Brown Brown, L.D. (1986) Fundamentals  of  statistical  exponential  families:  with  applications  in
statistical decision theory. Institute of Mathematical Statistics, Hayworth, CA, USA.
[Bühlman & van de Geer (2011)]BuhlmanDeGeer Bühlmann, P. & van de Geer, S. (2011). Statistics for High Dimensional Data: Methods, Theory and Applications. Springer Verlag.
[Candès & Plan (2009)]Candes Candès, E. & Plan, Y. (2009). Near-ideal model selection by ℓ_1 minimization. The Annals of Statistics, 37:5, 2145–2177.
[Clarke et al. (2009)]ClarkeEtal Clarke, B.S., Fokoué, E. & Zhang, H.H. (2009). Principles and Theory for Data Mining and Machine Learning. Springer Verlag.
[Cortes & Vapnik (1995)]Cortes Cortes, C. & Vapnik, V. (1995). Support-vector networks. Machine Learning 20 273–297.
[Cybenko (1989)]CybenkoCybenko, G. (1989). Approximation by Superpositions of a Sigmoidal Function 1989 Mathematics of Control,  Signals, and Systems, 2, 303–314.
[Darmois (1935)]Darmois Darmois, G. (1935). Sur les lois de probabilites a estimation exhaustive. Comptes Rendus de l'Académie des Sciencs, Paris, 200 1265–1266.
[Daubechies et al. (2004)]Daubechies  Daubechies, I. Defrise, M. & De Mol, C. (2004). An iterative thresholding algorithm for linear inverse problems with sparsity constraint. Communications on Pure and Applied Mathematics, 57:11, 1413–1457 
[Davison (1997)]Davison Davison, A.C. (1997). Bootstrap. Cambridge University Press.
[Davidson & MacKinnon (1993)]DavidsonMacKinnon1 Davidson, R. & MacKinnon, J.G. (1993). Estimation and Inference in Econometrics. Oxford University Press.
[Davidson & MacKinnon (2003)]DavidsonMacKinnon2 Davidson, R. & MacKinnon, J.G. (2003). Econometric Theory and Methods. Oxford University Press.
[Duo (1993)]Duo Duo, Q. (1993). The Formation of Econometrics. Oxford University Press.
[Debreu (1986)]Debreu Debreu, G. 1986. 
Theoretic Models: Mathematical Form and Economic Content. Econometrica, 54, 1259–1270.
[Dhillon et al. (2014)]Dhilonetal2013 Dhillon, P., Lu, Y. Foster, D.P. & Ungar, L.H. (2014). New Subsampling Algorithms for Fast Least Squares Regression. in Advances in Neural Information Processing Systems 26, Burges, Bottou, Welling, Ghahramani & Weinberger Eds., Curran Associates.
[Efron & Tibshirani (1993)]Efron Efron, B. & Tibshirani, R. (1993). Bootstrap. Chapman Hall CRC.
[Engel (1857)]Engel Engel, E. (1857). Die Productions- und Consumtionsverhältnisse des Königreichs Sachsen. Statistisches Bureau des Königlich Sächsischen Ministeriums des Innern.
[Feldstein & Horioka (1980)]FH Feldstein, M. & Horioka, C. (1980). Domestic Saving and International Capital Flows. Economic Journal, 90, 314–329.
[Flach (2012)]Flach Flach, P. (2012). Machine Learning. Cambridge University Press.
[Foster & George (1994)]FosterGeorge94 Foster, D.P. & George, E.I. (1994). The Risk Inflation Criterion for Multiple Regression. The Annals of Statistics, 22:4, 1947–1975.
[Friedman (1997)]Friedman97 Friedman, J.H. (1997). Data Mining and Statistics: What's the Connection. Proceedings of the 29th Symposium on the Interface Between Computer Science and Statistics.
[Frish & Waugh (1933)]FrishWaugh Frisch, R. & Waugh, F.V. (1933). Partial Time Regressions as Compared with Individual Trends. Econometrica. 1, 387–401.
[Gneiting (2011)]Gneiting Gneiting, T. (2011). Making and Evaluating Point Forecasts. Journal of the American Statistical Association, 106, 746–762.
[Givord (2010)]Givord Givord, P. (2010). Méthodes économétriques pour l'évaluation de politiques publiques. INSEE Document de Travail, 08
[Grandvalet et al. (2005)]Grandvalet Grandvalet, Y., Mariéthoz, J., & Bengio, S. 2005. Interpretation of SVMs with an application to unbalanced classification. Advances in Neural Information Processing Systems 18.
[Groves & Rothenberg (1969)]GrovesRothenberg69 Groves, T. & Rothenberg, T. (1969). A note on the expected value of an inverse matrix. Biometrika, 56:3, 690–691.
[Haavelmo (1944)]Haavelmo Haavelmo, T. (1944). The probability approach in econometrics, Econometrica, 12:iii-vi and 1–115.
[Hastie & Tibshirani (1990)]Trevor Hastie, T. & Tibshirani, R. (1990). Generalized Additive Models. Chapman & Hall/CRC.
[Hastie et al. (2009)]HastieEtal Hastie, T., Tibshirani, R. & Friedman, J. (2009). The Elements of Statistical Learning. Springer Verlag.
[Hastie (2005)]Hastie Hastie, T., Tibshirani, W. & Wainwright, M. (2015). Statistical Learning with Sparsity. Chapman CRC.
[Hastie et al. (2016)]HTT  Hastie, T., Tibshiriani, R. & Tibshiriani, R.J. (2016). Extended comparisons of best subset selection, forward stepwise selection and the Lasso. ArXiV, https://arxiv.org/abs/1707.08692.
[d'Haultefœuille & Givord (2014)]dHaultefoeuille  d'Haultefœuille, X. & Givord, P. (2014) La régression quantile en pratique. Économie & Statistiques, 471, 85–111.
[Hebb (1949)]Hebb Hebb, D.O. (1949). The organization of behavior, New York, Wiley.
[Heckman (1979)]Heckman Heckman, J.J. (1979). Sample selection bias as a specification error. Econometrica, 47, 153–161.
[Heckman et al. (2003)]Heckman2 Heckman, J.J., Tobias, J.L. & Vytlacil, E. (2003). Simple Estimators for Treatment Parameters in a Latent-Variable Framework. The Review of Economics and Statistics, 85, 748–755.

[Hendry & Krolzig (1995)]HendryK Hendry, D F. & Krolzig, H.-M. (2001). Automatic Econometric Model Selection. Timberlake Press.
[Herbrich et al. (1999)]Herbrich Herbrich, R., Keilbach, M., Graepel, T. Bollmann-Sdorra, P. & Obermayer, K. (1999). Neural Networks in Economics. in Computational Techniques for Modelling in Economics, T. Brenner Eds. Springer Verlag, 169–196.
[Hoerl (1962)]Hoerl Hoerl, A.E. (1962). Applications of ridge analysis to regression problems. Chemical Engineering Progress, 58:3, 54–59.
[Hoerl & Kennard (1980)]HoerlKennard Hoerl, A.E. & Kennard, R.W. (1981). Ridge regression: biased estimation for nonorthogonal problems This Week's Citation Classic, ISI, http://bit.ly/2H9LGiD
[Holland (1986)]Holland Holland, P. (1986). Statistics and causal inference. Journal of the American Statistical Association, 81, 945–960.
[Hyndman et al. (2009)]Hyndman Hyndman, R. , Koehler, A.B., Ord, J.K.	& Snyder, R.D. (2009). Forecasting with Exponential Smoothing. Springer Verlag.
[James et al. (2013)]James James, G., D. Witten, T. Hastie, & R. Tibshirani (2013). An introduction
to Statistical Learning. Springer Series in Statistics.
[Khashman (2011)]Khashman Khashman, A. (2011). Credit risk evaluation using neural networks: Emotional versus conven-
tional models. Applied Soft Computing, 11, 5477–5484.
[Kean (2010)]Keen Kean, M.P. (2010). Structural vs. atheoretic approaches to econometrics. Journal of Econometrics, 156, 3–20.
[Kleiner et al. (2012)]  Kleiner, A., Talwalkar, A., Sarkar , P. & Jordan, M. (2012). The Big Data Bootstrap. arXiv:1206.6415 .
[Koch (2013)]Koch Koch, I. (2013). Analysis of Multivariate and High-Dimensional Data. Cambridge University Press. 
[Koenker (1998)]KoenkerGalton Koenker, R. (1998). Galton, Edgeworth, Frish, and prospects for quantile regression in Econometrics. Conference on Principles of Econometrics, Madison.
[Koenker (2003)]Koenker Koenker, R. (2003). Quantile Regression. Cambridge University Press.
[Koenker & Machado (1999)]KoenkerMachado Koenker, R. & Machado, J. (1999). Goodness of fit and related inference processes for quantile regression Journal of the American Statistical Association, 94, 1296-1309.
[Kolda & Bader (2009)]Kolda Kolda, T. G. & Bader, B. W. (2009). Tensor decompositions and applications. SIAM Review 51, 455–500.
[Koopmans (1957)]Koopmans Koopmans, T.C. (1957). Three   Essays   on   the   State   of   Economic   Science. McGraw-Hill.
[Kuhn & Johnson (2013)]KuhnJohnson Kuhn, M. & Johnson, K. (2013). Applied Predictive Modeling. Springer Verlag.
[Landis & Koch (1977)]LandisKoch Landis, J.R. & Koch, G.G. (1977). The measurement of observer agreement for categorical data. Biometrics, 33, 159–174. 
[LeCun et al. (2015)]LeCun LeCun, Y., Bengio, Y. & Hinton, G. (2015). Deep learning. Nature 521 436–444.
[Leeb (2008)]Leeb Leeb, H. (2008). Evaluation and selection of models for out-of-sample prediction when the sample size is small relative to the complexity of the data-generating process. Bernoulli 14:3, 661–690.
[Lemieux (2006)]Lemieux Lemieux, T. (2006). The « Mincer Equation » Thirty Years After Schooling, Experience, and Earnings. in  Jacob Mincer A Pioneer of Modern Labor Economics, Grossbard Eds, 127–145, Springer Verlag.
[Li & Racine (2006)]Li Li, J. & J. S. Racine (2006). Nonparametric Econometrics. Princeton University Press.
[Li et al. (2017)]Lili Li, C., Li, Q., Racine, J. & Zhang, D. (2017). Optimal Model Averaging Of Varying Coefficient Models. Department of Economics Working Papers 2017-01, McMaster University.
[Lin et al. (2016)]Lin Lin, H.W., Tegmark, M. & Rolnick, D. (2016). Why does deep and cheap learning work so well?
ArXiv e-prints.
[Lucas (1976)]Lucas Lucas, R.E. (1976). Econometric Policy Evaluation: A Critique. Carnegie-Rochester Conference Series on Public Policy, 19–46.
[Mallows (1973)]Mallows Mallows, C.L. (1973). Some Comments on C_p. Technometrics, 15, 661–675.
[McCullogh & Pitts (1943)]McCulloch McCullogh, W.S. & Pitts, W. (1943). A logical calculus of the ideas immanent in nervous activity. Bulletin of Mathematical Biophysics, 5:4, 115–133.
[Mincer (1974)]Mincer Mincer, J. (1974). Schooling, experience and earnings. Columbia University Press.
[Mitchell (1997)]Mitchell Mitchell, T. (1997). Machine Learning. McGraw-Hill.
[Morgan & Sonquist (1963)]Morganetal Morgan, J.N. & Sonquist, J.A. (1963). Problems in the analysis of survey data, and a proposal. Journal of the American Statistcal Association, 58, 415–434.
[Morgan (1990)]Morgan Morgan, M.S. (1990). The history of
econometric ideas. Cambridge University Press.
[Mohri et al. (2012)]MohriEtal Mohri, M., Rostamizadeh, A. & Talwalker, A. (2012) Foundations of Machine Learning. MIT Press.
[Mullainathan & Spiess (2017)]Mullainathan
Mullainathan, S. & Spiess, J. (2017). Machine learning: An applied econometric approach. Journal of Economic Perspectives, 31 87–106.
[Müller (2011)]Muller Müller, M. (2011). Generalized Linear Models in  Handbook of Computational Statistics, J.E Gentle, W.K. Härdle & Y. Mori Eds. Springer Verlag.
[Murphy (2012)]Murphy Murphy, K.R. (2012). Machine Learning: a Probabilistic Perspective. MIT Press.
[Murphy & Welch (1990)]MurphyWelch Murphy, K. M. & F. Welch (1990). Empirical age-earnings profiles. Journal of Labor Economics 8, 202–229.
[Nadaraya (1964)]Nadaraya Nadaraya, E. A. (1964). On Estimating Regression. Theory of Probability and its Applications, 9:1, 141–2.
[Natarajan (1995)]Natarajan  Natarajan, B. K. (1995). Sparse approximate solutions to linear systems. SIAM Journal on Computing (SICOMP), 24 227–-234.
[Nevo & Whinston (2010)]Nevo Nevo, A. & Whinston, M.D. (2010). Taking the Dogma out of Econometrics: Structural Modeling and Credible Inference. Journal of Economic Perspective, 24,  69–82.
[Neyman (1923)]Neyman Neyman, J. (1923).  Sur les applications de la théorie des probabilités aux expériences agricoles : Essai des principes. Mémoire de master, republibé dans Statistical Science, 5, 463–472.
[Nisbet, Elder & Miner (2001)]Nisbet Nisbet, R., Elder, J. & Miner, G. (2011). Handbook of Statistical Analysis and Data Mining
Applications. Academic Press, New York. 
[Okun (1962)]Okun Okun,  A. (1962). Potential GNP: Its measurement and significance. Proceedings of the Business and Economics Section of the American Statistical Association, 98–103.
[Orcutt (1952)]Orcutt Orcutt, G.H. (1952). Toward a partial redirection of econometrics. Review of Economics and Statistics, 34 195–213.
[Pagan & Ullah (1999)]Pagan Pagan, A. & A. Ullah (1999). Nonparametric Econometrics. Themes in Modern Econometrics. Cambridge: Cambridge University Press.
[Pearson (1901)]Pearson Pearson, K. (1901). On lines and planes of closest fit to systems of points in space. Philosophical Magazine, 2, 559-–572.
[Platt (1999)]Platt Platt, J. (1999). Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. Advances in Large Margin Classifiers. 10, 61–74.
[Portnoy (1988)]Portnoy Portnoy, S. (1988). Asymptotic Behavior of Likelihood Methods for Exponential Families when the Number of Parameters Tends to Infinity. Annals of Statistics, 16:356–366.
[Quenouille (1949)]Quenouille1 Quenouille, M. H. (1949). Problems in Plane Sampling. The Annals of Mathematical Statistics 20(3):355–375.
[Quenouille (1956)]Quenouille2 Quenouille, M. H. (1956). Notes on Bias in Estimation. Biometrika 43(3-4), 353–360.
[Quinlan (1986)]Quinlan Quinlan, J.R. (1986). Induction of decision trees. Machine Learning  1 81–106.
[Reiersøol (1945)]Reiersol Reiersøol, O. (1945). Confluence analysis of means of instrumental sets of variables. Arkiv. for 
Mathematik, Astronomi Och Fysik, 32. 
[Rosenbaum & Rubin (1983)]RosenbaumRubin Rosenbaum, P. & Rubin, D. (1983). The Central Role of the Propensity Score in Observational Studies for Causal Effects. Biometrika, 70, 41–55.
[Rosenblatt (1958)]Rosenblatt Rosenblatt, F. (1958). The perceptron: a probabilistic model for information storage and organization in the brain. Psychological Review, 65, 386–408.
[Rubin (1974)]Rubin Rubin, D. (1974). Estimating Causal Effects of Treatments in Randomized and Nonrandomized Studies. Journal of Educational Psychology, 66, 688–701.
[Ruppert, Wand & Carroll (2003)]Ruppert Ruppert, D., Wand, M. P. & Carroll, R.J. (2003). Semiparametric Regression. Cambridge University Press.
[Samuel (1959)]Samuel Samuel, A. (1959). Some Studies in Machine Learning Using the Game of Checkers. IBM Journal of Research and Development, 44:1. 
[Schultz (1930)]Schultz Schultz, H. (1930). The Meaning of Statistical Demand Curves. University of Chicago.
[Shai & Shai (2014)]Shaix2 Shai, S.S. & Shai, B.D. (2014). Understanding Machine Learning
From Theory to Algorithms. Cambridge University Press.
[Shao (1993)]Shao93 Shao, J. (1993). Linear Model Selection by Cross-Validation. Journal of the American Statistical Association 88:(422), 486–494.
[Shalev-Shwartz &  Ben-David (2014)]Shalev-Shwartz Shalev-Shwartz, S. &  Ben-David, S. (2014). Understanding Machine Learning: From Theory to Algorithms. Cambridge University Press.
[Shao (1997)]Shao97 Shao, J. (1997). An Asymptotic Theory for Linear Model Selection. Statistica Sinica, 7, 221–264.
[Shapire & Freund (2012)]ShapireFreund Shapire, R.E. & Freund, Y. (2012). Boosting. MIT Press.
[Silverman (1986)]Silverman Silverman, B.W. (1986) Density Estimation. Chapman & Hall.
[Simonoff (1996)]Simonoff Simonoff, J. S. (1996). Smoothing Methods in Statistics. Springer. 
[Stone (1977)]Stone Stone, M. (1977). An Asymptotic Equivalence of Choice of Model by Cross-Validation and Akaike's Criterion. Journal of the Royal Statistical Society. Series B , 39:1, 44–47.
[Tam & Kiang (1992)]Tam Tam, K.Y. & Kiang, M.Y. (1992). Managerial applications of neural networks: The case
of bank failure predictions. Management Science, 38, 926–947.
[Tan (1995)]Tan Tan, H. (1995). Neural-Network model for stock forecasting. MSc Thesis, Texas Tech. University.
[Tibshirani (1996)]LASSO Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society, Series B., 58, 267–288.
[Tibshirani & Wasserman (2016)]TibWasserman Tibshirani, R. & Wasserman, L. (2016). A Closer Look at Sparse Regression. http://bit.ly/2FrGQ32
[Tikhonov (1963)]Tikhonov Tikhonov, A. N. (1963). Solution of incorrectly formulated problems and the regularization method. Soviet Mathematics, 4: 1035–1038.
[Tinbergen (1939)]Tinbergen Tinbergen, J. (1939). Statistical Testing of Business Cycle Theories. Vol. 1: A Method and its Application to Investment activity; Vol. 2: Business Cycles in the United States of America,
1919—1932. Geneva: League of Nations.
[Tobin (1958)]Tobin Tobin, J. (1958). Estimation of Relationship for Limited Dependent Variables. Econometrica, 26, 24–36. 
[Tropp (2011)]Tropp2011 Tropp, (2011). Improved analysis of the subsampled randomized Hadamard transform. Advances in Adaptive Data Analysis, 3:1, 115–126.
[Tsen (2001)]Tsen Tsen, P. (2001). Convergence of a block coordinate descent for nondifferentiable minization. Journal of Optimization Theory and Applications, 109:3, 475–494.
[Tufféry (2001)]Tuffery Tufféry, S. (2001). Data Mining and Statistics for Decision Making. Wiley Interscience.
[Tukey (1958)]Tukey Tukey, J. W. (1958). Bias and confidence in not quite large samples. The Annals of Mathematical Statistics, 29:614–623.
[Vapnik (1998)]Vapnik98 Vapnik, V. (1998). Statistical Learning Theory. Wiley.
[Vapnik & Chervonenkis (1971)]VC Vapnik, C, & Chervonenkis, A. (1971). On the uniform convergence of relative frequencies of events to their probabilities. Theory of Probability and its Applications, 16:264–280.
[Varian (2014)]Varian Varian, H.R. (2014). Big Data: New Tricks for Econometrics. Journal of Economic Perspectives,  28(2):3–28. 
[Vert (2017)]Vert Vert, J.P. (2017). Machine learning in computational biology. ENSAE.
[Waltrup et al. (2014)]Walltrup Waltrup, L.S., Sobotka, F., Kneib, T. & Kauermann, G. (2014). Expectile and quantile regression—David and Goliath? Statistical Modelling, 15, 433 – 456.
[Watson (1964)]Watson Watson, G. S. (1964). Smooth regression analysis. Sankhya: The Indian Journal of Statistics, Series A, 26:4, 359–372.
[Watt et al. (2016)]Watt Watt, J., Borhani, R. & Katsaggelos, A. (2016). Machine Learning Refined : Foundations, Algorithms, and Applications. Cambridge University Press.
[Widrow & Hoff (1960)]Widrow Widrow, B. & Hoff,  M.E. Jr. (1960). Adaptive Switching Circuits. IRE WESCON Convention Record, 4:96–104.
[Wolpert & Macready (1997)]Wolpert97 Wolpert, D.H., Macready, W.G. (1997), No Free Lunch Theorems for Optimization, IEEE Transactions on Evolutionary Computation 1, 67.
[Wolpert (1996)]Wolpert96 Wolpert, David (1996), The Lack of A Priori Distinctions between Learning Algorithms, Neural Computation, 1341-1390.
[Working (1927)]Working Working, E. J. (1927). What do statistical `demand curves' show? Quarterly Journal of Economics, 41:212–35.
[Yu & Moyeed (2001)]YuMoyeed Yu, K. & Moyeed, R. (2001). Bayesian quantile regression. Statistics & Probability Letters, 54, 437–447.

[Zinkevich et al. (2010)]Zinkevichetal Zinkevich M.A., Weimer, M., Smola, A. & Li, L. (2010). Parallelized Stochastic Gradient. Advances in neural information processing systems, 2595–2603.


