# Concilier risques collectifs et décisions individuelles

cet article a été co-écrit avec Laurence Barry.

Les débuts de la pandémie de SARS-CoV-2 (ou COVID-19) ont vu se multiplier les appels à la « responsabilité individuelle », en commençant par de fortes demandes (voire une obligation dans certains pays, dont la France) à rester chez soi autant que possible, au début du printemps 2020, avant qu’il ne soit obligatoire de porter un masque dans les lieux publics (souvent fermés) au cours de l’été. En paraphrasant Coluche, « dire qu’il suffirait que les gens restent chez eux pour qu’on puisse sortir… ». Cet appel à la responsabilité de chacun est faite au nom de tous et pour le bien de tous, venant symboliser cette solidarité toute particulière que nous rappelle la pandémie : le risque que je choisis de courir ne concerne pas seulement ma personne mais constitue aussi un risque pour ceux qui m’entourent. Pour le formuler en terme probabiliste, McKendrick (1926) affirmait « la probabilité d’occurrence augmente avec le nombre de cas existants ». Assez intuitive a priori, cette conception de la responsabilité individuelle va en réalité à l’encontre de la conception classique en économie : l’individu rationnel (et responsable) fait des choix le concernant, et ne concernant que lui. Le bien collectif se déduit par sommation des utilités individuelles, indépendantes les unes des autres. Seulement voilà ; avec l’épidémie se crée une interdépendance des utilités qui fait que le bien-être d’Untel, qui choisit de ne pas porter de masque, peut nuire à la santé et donc l’utilité de beaucoup d’autres personnes. Comment penser alors en termes économiques cette « responsabilité individuelle » dans le contexte de l’épidémie ?

# Des préférences individuelles au bien-être collectif

L’hypothèse centrale de la théorie économique du comportement est que chacun est capable de classer, par ordre de préférence, toutes sortes d’alternatives qui lui sont proposées. Et si je dois choisir une parmi deux, je choisirai systématiquement celle que je préfère. Comme le montre Mas-Colell et al. (1995), une simple hypothèse de continuité des préférences se traduit alors par l’existence d’une fonction d’utilité individuelles traduisant ces préférences. Cette approche pourrait suffire dans l’état de nature de Jean-Jacques Rousseau, lorsque l’homme est imaginé vivant en solitaire. Mais en société il convient d’être plus réaliste, et de tenir compte des interactions entre les individus. Organiser la vie en société, en favorisant la coopération et en cherchant à assurer un bien-être collectif, ne peut en effet se faire en se contentant de comprendre le bien-être individuel. Pour reprendre un exemple de Jean-Jacques Rousseau, plusieurs chasseurs ont intérêt à collaborer pour traquer un cerf, car aucun chasseur ne saurait y arriver seul[i] (Rousseau, 1755). La première difficulté est donc d’assurer une collaboration pour la chasse, mais aussi et surtout, si un cerf est tué, se pose le problème de la répartition de la viande.

Tout au long du XVIIIème siècle Francis Hutcheson et Adam Smith en Angleterre, Jean-Charles de Borda et Nicolas de Condorcet ont tenté de formaliser cette notion de « bien-être collectif », montrant qu’il existait malheureusement de très nombreux paradoxes, en particulier quand il s’agit du bien être d’une nation. A la même époque, Kant formalise l’impératif catégorique : « Agis uniquement d’après la maxime qui fait que tu puisses vouloir en même temps qu’elle devienne une loi universelle ». Autrement dit, avant de prendre une décision pour agir, il convient de se demander ce qui se passerait si tout le monde agissait de cette manière. La rationalité de l’individu se doit d’être collective, et de prendre en compte l’humanité dans son ensemble. Dans sa lecture de Kant, Arendt (1991) met en avant le sensus communis, ce sens commun à tous les hommes et qui rattache le jugement de chacun « à la raison humaine tout entière ». Penser par soi-même devient alors « penser en se mettant à la place de tout autre, dans ce qu’elle appelle une « mentalité élargie ».

Tocqueville quant à lui renverse les termes de l’équation. Dans les pays démocratiques selon lui, on ne peut plus mettre en avant la valeur du sacrifice de soi : il faut pouvoir démontrer que « l’homme en servant ses semblables se sert lui-même (…) Aux Etats-Unis on ne dit presque point que la vertu est belle. On dit qu’elle est utile ». Cela implique cependant de « petits sacrifices », consentis car ils se révèlent bénéfiques pour celui qui les consent. Tocqueville (1981) exhorte alors ses lecteurs à agir dans leur intérêt « bien entendu », c’est-à-dire en tenant compte de l’intérêt de tous.

Cette rationalité collective est en réalité une pratique habituelle au sein de petits groupes, comme la famille. Il n’est pas rare, en effet, de mettre de côté son intérêt personnel pour le bien de la famille. Mais elle est plus complexe à mettre en œuvre au sein d’un groupe plus important, plus hétérogène, voire plus abstrait.

# Le passager clandestin contre les intérêts communs

Une vision résolument optimiste consisterait en effet à croire que si tous les membres d’un groupe ont des intérêts communs, alors chacun va agir pour les atteindre. Un exemple bien connu est celui du réchauffement climatique : collectivement, l’intérêt de tous est la réduction des gaz à effet de serre au niveau mondial ; mais individuellement, chaque pays a la tentation de retarder la mise en place de mesures qui pourraient pénaliser son économie, en espérant toutefois bénéficier d’actions précoces de pays voisins. C’est le principe du passager clandestin : il y aurait un bénéfice collectif à tirer d’une coopération, mais les individus ont davantage d’incitations à chercher à profiter de la « coopération » des autres. En termes économiques, ils cherchent à avoir une prestation sans en assumer les coûts.

Ce problème, largement étudié dans Olson (1965), est classique pour la majorité des « biens publics[ii] » qui satisfont deux caractéristiques : être non-rival et non-excluable, c’est-à-dire dont la consommation par les uns ne diminue pas la quantité disponible pour les autres et dont on ne peut par ailleurs restreindre l’accès. Axelrod & Hamilton (1981) expliquaient que la coopération nécessaire à la promotion de biens communs ne dépendait pas forcément d’une forme d’altruisme, mais plus simplement d’une réciprocité entre les agents, basée sur une coopération conditionnelle : ils coopèrent s’ils pensent que les autres vont faire de même. Plusieurs études ont montré que la majorité des gens fonctionnent de la sorte, mais leur comportement est très sensible à leurs croyances, d’où l’importance de maintenir leurs convictions en matière d’égalité (ou d’égalitarisme) : tout le monde doit coopérer, personne ne doit bénéficier d’un traitement de faveur. Fehr & Fischbacher (2004) ont ainsi montré qu’il suffit d’une petite proportion de passagers clandestins pour provoquer une rupture[iii] de la coopération. Ceci explique probablement les diverses normes injonctives autour de la « distanciation sociale », assurant qu’une personne qui resquille soit sanctionnée de manière exemplaire. En effet, comme le soulignent Brito et al. (1991) à propos des vaccins, même si l’obligation de vacciner est sous-optimale, elle peut être nécessaire si la proportion de gens prêts à se porter volontaires est en dessous du seuil nécessaire à l’immunité de la population dans son ensemble.

# Le cas de la vaccination

La vaccination est en fait un exemple presque parfait de ce problème de passager clandestin, via la notion d’immunité de groupe. Plus le taux de personnes immunisées augmente, dans un groupe, plus le risque pour une personne non-immunisée de rencontrer une personne infectieuse diminue, et au-delà d’un certain seuil (de l’ordre de 80% pour la plupart des maladies, comme la coqueluche, la variole, la polio, etc), il devient impossible pour la maladie de se maintenir dans la population et elle finit par disparaître. Pour les maladies contagieuses bénéficiant d’un vaccin, au niveau collectif, il est souhaitable que 80% de la population soit vaccinée ; mais si la vaccination a des effets secondaires conséquents, il peut être rationnel au niveau individuel de ne pas souhaiter être vacciné.

Deux aspects importants entrent alors en jeu : la perception des risques, et la croyance dans le comportement des risques des autres membres de la communauté. Avoir une minorité de passagers clandestins (disons moins de 20%), parce qu’ils pensent les risques trop grands, n’est pas problématique. Mais si la perception des risques change, on peut observer la rupture de l’équilibre, et l’immunité de groupe n’existe plus. Aussi, la confiance dans l’autorité est essentielle, comme le rappelait Charpentier (2020).

L’immunité collective fonctionne grâce à un contrat social implicite : ceux qui sont médicalement capables de se faire vacciner doivent se faire vacciner. La contrepartie est que les personnes que ne souhaitent pas respecter ce contrat devraient s’engager à ce que leurs actions n’entraînent pas de coût supplémentaire pour ceux qui le respectent, en particulier en s’imposant une forte distanciation sociale, en évitant les lieux publics, de manière à ne pas contaminer des personnes ayant de faibles défenses immunitaires, et qui sont, elles, dans l’obligation de compter sur l’immunité collective.

# Les pandémies et leurs réponses individuelles

Comme le disait Daniel Kahneman dans Konnikova (2020), « people, certainly including myself, don’t seem to be able to think straight about exponential growth. What we see today are infections that occurred 2 or 3 weeks ago and the deaths today are people who got infected 4 or 5 weeks ago. All of this is I think beyond intuitive human comprehension ». Le fait d’adopter une attitude de passager clandestin et de ne pas respecter les contraintes de distanciation sociale tient peut-être simplement du fait qu’on ne comprend simplement pas ce qu’est une croissance exponentielle : on ne mesure pas vraiment l’impact de sa propre contagion sur le groupe dans son ensemble.

Ce point a été montré dans Lammers et al. (2020) à partir de l’interprétation du nombre de reproduction de base R0 des modèles épidémiologiques. Le  Rcorrespond au nombre moyen de personnes qu’une personne contagieuse peut infecter : il constitue ainsi une visualisation de la contagiosité d’un individu sur son entourage. Avec un R0 de 1, la pandémie est contrôlée, et la croissance est linéaire. Mais s’il excède 1, la croissance est exponentielle. Avec un R0 de 1.5, 4 personnes vont en contaminer 6 autres, qui à leur tour vont en contaminer 9 autres, etc. En une quinzaine d’itérations, 1 750 personnes seront contaminées. Avec un R0 de 2, ces 4 individus auront contaminé plus de 130 000 personnes (soit 75 fois plus), en une quinzaine d’itération ! Autrement dit, alors qu’au niveau de l’individu spécifique l’augmentation est à peine perceptible (il contamine 2 personnes au lieu de 1.5 en moyenne), l’effet collectif est, lui, extrêmement important et difficilement concevable.

Figure 1 : nombre de personnes contaminées après 1, 3, 5, 7, 9 itérations, pour différentes valeurs de R0 (entre 1.8 en haut à gauche et 2.4 en bas à droite).

De nombreuses études en science du comportement ont montré que nous sommes davantage sensibilisés face à une seule personne identifiable qu’en étant noyé sur une avalanche de chiffres. Ce serait alors par l’exemple que l’on pourrait se convaincre mutuellement de coopérer pour le bien commun. Le port du masque facial est intéressant, car si des sondages ont montré qu’une majorité des gens portaient un masque pour se protéger, les masques ont surtout pour effet de protéger les autres personnes d’une transmission asymptomatique du SARS-CoV-2. Il présente aussi l’avantage de rendre visible la nouvelle norme sociale et d’impliquer activement tous les membres de la communauté. En devenant un symbole de la solidarité, le port du masque engage ainsi la coopération de chacun. A l’inverse, les photos de personnes à la plage ou dans les parcs publics qui ne respectent pas la distanciation sociale ont probablement eu un impact préjudiciable en termes de changement de comportement.

Les débats autour des applications de traçage sont un autre exemple frappant de la difficulté de faire accepter la coopération. Comme pour le masque, elles sont présentées le plus souvent comme permettant d’être alerté si l’on a été en contact avec une personne contaminée, donc comme un moyen de se protéger soi-même ; beaucoup plus rarement est mise avant la possibilité de prévenir autrui de sa propre contamination, parfois un inconnu qu’on ne pourra jamais alerter sans l’application. De plus, alors qu’elles étaient recommandées par de nombreux épidémiologistes (di Domenico et al. 2020; Ferguson et al. 2020; Ferretti et al. 2020), elles ont été dénoncées soit parce qu’elles porteraient atteinte à la liberté individuelle, soit parce qu’elles présenteraient des dangers de détournement. Dans une importante contribution, des experts en cryptographie ont ainsi tenté d’alerter l’opinion publique sur les possibles usages malveillants de ces applications ; au travers d’une quinzaine d’exemples qui cherchent à marquer l’imagination, ce n’est plus le passager clandestin qui est mis en avant pour saper la coopération mais l’individu franchement malveillant qui chercherait à nuire à ses voisins (Vuillot et al. 2020).

# De la place de la liberté individuelle

Dans le contexte de la vaccination, et plus récemment sur le port du masque, l’argument de la liberté de choix est souvent avancé par les opposants, laissant croire que l’exercice de la liberté se faisait sans contraintes. Comme le note Frankfurt (2003), la plupart des religions limite le comportement d’une personne dans la mesure où elle agit en accord avec les préceptes de son Dieu ou de son église. Dans un contexte plus laïque, Jean-Jacques Rousseau affirmait que « l’obéissance à la loi qu’on s’est prescrite est liberté ». Car le concept de liberté s’accompagne toujours de la notion de responsabilité : je suis libre lorsque deux conditions sont réunies. Premièrement, j’ai la capacité d’agir (ou de ne pas agir) d’une manière particulière, et deuxièmement, j’accepte la responsabilité de mes actes. Quand je refuse de porter un masque en période de pandémie, j’accepte la première condition mais rejette la seconde. Autrement dit, j’affirme mon droit d’agir, ou de ne pas agir, mais je le fais de manière à refuser d’accepter toute responsabilité pour les conséquences que mes actions (ou inactions) peuvent entraîner. Comme l’affirmait Friedrich Hayek, « la liberté ne signifie pas seulement qu’une personne a le droit de choisir et qu’elle porte le fardeau de ses choix, mais aussi qu’elle doit assumer les conséquences de ses actes, pour lesquels elle sera félicitée ou blâmée. Liberté et responsabilité sont indissociables ».

# L’assurance comme réponse collective ?

Dans le contexte de la santé, Wikler (2002) affirmait « if people know they are taking risks but accept them as the price of pursuing goals to which they assign higher priority, then it is not the business of public health to insist that health be valued above all ». Si ce précepte peut valoir pour l’assurance santé classique, il est plus difficile à appliquer à l’épidémie ; comme expliqué plus haut, dans les maladies contagieuses le choix individuel de prendre un risque se répercute sur le reste de la collectivité. De plus, la logique assurantielle de couverture de l’aléa grâce à la mutualisation fonctionne mal dans le cadre de l’épidémie ; on est dans le cas classique d’un risque systémique où les individus et leurs risques ne sont pas indépendants. En réalité, la contagion met en avant une solidarité d’un autre ordre que celle promue par l’assurance, comme le note Barry (2020). Il s’agit d’une interdépendance où le comportement de l’un impacte le risque de l’autre et qui exige, pour être contrôlé, la coopération de tous. Penser collectivement, c’est donc finalement adopter des valeurs de solidarité et de coopération.

# Références

Arendt, Hannah (1991). Juger – La Philosophie Politique de Kant. Points. Paris: Seuil.

Axelrod, Robert & Hamilton, William (1981). The Evolution of Cooperation. Science, 211 (4489): 1390–96.

Barry, Laurence (2020). Individu/Collectif : L’épidémiologie à l’épreuve Du Big Data (Ou l’inverse) ? Working Paper # 20. Paris: Chaire PARI.

Brito, Dagobert, Sheshinski Eytan & Intriligator Michael (1991). Externalities and Compulsory Vaccinations. Journal of Public Economics 45, 69–90.

Charpentier, Arthur (2020). De la démarche scientifique en période de crise. Risques, 121.

Costa, Dora et Kahn, Matthew (2003). Civic Engagement and Community Heterogeneity: An Economists Perspective. Perspectives on Politics, 1(1): 103-112.

di Domenico, Laura, Giulia Pullano, Chiara Sabbatini, Pierre-Yves Boelle, and Vittoria Colizza. (2020). Expected Impact of Lockdown in Île-de-France and Possible Exit Strategies. 9. Paris: INSERM.

Fehr, Ernst & Fischbacher, Urs (2004). Social norms and human cooperation. Trends in Cognitive Sciences, 8(4), 185-190.

Ferguson, Neil, D. Laydon, G. Nedjati Gilani, N. Imai, K. Ainslie, M. Baguelin, S. Bhatia, A. Boonyasiri, Z. Cucunuba Perez, G. Cuomo-Dannenburg, A. Dighe, I. Dorigatti, H. Fu, K. Gaythorpe, W. Green, A. Hamlet, W. Hinsley, L. Okell, S. Van Elsland, H. Thompson, R. Verity, E. Volz, H. Wang, Y. Wang, P. Walker, C. Walters, P. Winskill, C. Whittaker, C. Donnelly, S. Riley, & A. Ghani. 2020. Report 9: Impact of Non-Pharmaceutical Interventions (NPIs) to Reduce COVID19 Mortality and Healthcare Demand. Imperial College Report.

Ferretti, Luca, Chris Wymant, Michelle Kendall, Lele Zhao, Anel Nurtay, Lucie Abeler-Dörner, Michael Parker, David Bonsall, & Christophe Fraser. (2020). Quantifying SARS-CoV-2 Transmission Suggests Epidemic Control with Digital Contact Tracing. Science, Vol. 68 #619, 1-8.

Frankfurt, Harry, (2003) Freedom of the Will and a Concept of a Person, in Gary Watson (ed), Free Will, 2nd edition, Oxford University Press,  322-336.

Konnikova, Maria (2020). Why We Underestimated COVID-19. New Yorker, 3 avril 2020,

Lammers, Joris, Crusius, Jan et Gast, Anne (2020). Correcting misperceptions of exponential coronavirus growth increases support for social distancing. PNAS, 117 (28).

Lim, Wooyoung & Zhang, Pengfei (2020). Herd immunity and a vaccination game: An experimental study. PLoS One. 5(5), e0232652.

McKendrick, A. G. (1926). Applications of Mathematics to Medical Problems. Proceedings of the Edinburgh Mathematical Society 44, 98–130.

Mas-Colell, Andreu, Whinston, Michael et Green, Jerry (1995). Microeconomic Theory. Oxford University Press.

Olson, Mancur (1965). The Logic of Collective Action: Public Goods and the Theory of Groups. Harvard University Press.

Rousseau, Jean-Jacques (1755). Discours sur l’origine et les fondements de l’inégalité parmi les hommes. Garnier Flammarion.

Skyrms, Brian (2004) The Stag Hunt and the Evolution of Social Structure. Cambridge: Cambridge University Press.

Tocqueville, Alexis. (1981). De La Démocratie En Amérique – 2. Flammarion.

Vuillot, Xavier, Anne Bonnetain, Veronique Canteaut, Pierrick Cortier, Lucca Gaudry, Steve Hirschi, Stéphanie Kremer, Gaëtan Lacour, Matthieu Leurent, Léo Lequesne, André Perrin, Emmanuel Schrottenloher, Serge Thomé, and Christophe Vaudenay. (2020). Le Traçage Anonyme, Dangereux Oxymore Analyse de Risques à Destination Des Non-Spécialistes. Https://Risques-Tracage.Fr. Retrieved April 25, 2020 (https://risques-tracage.fr/).

Wikler, Daniel (2002) Personal and social responsibility for health. Ethics & International Affairs, 16, 47-55.

Wolman, Davis (2020). Yes, the Public Can Be Trusted in a Pandemic. Wired, 27 mars 2020,

[i] Comme le rappelle Skyrms (2004), ce dilemme de chasse au cerf est aussi appelé jeu de l’assurance, en théorie des jeux et de choix sociaux.

[ii] On pourrait aussi penser aux « biens communs », qui sont des biens en capacité limité. Le SARS-CoV-2 a montré que le système de santé pouvait être saturé, ce qui en fait dès lors un bien rival.

[iii] Mathématiquement, cette rupture est intéressante car on peut alors voir la vaccination comme un jeu non-linéaire de bien public, comme le fait Lim & Zhang (2020).

# Gini index, poverty and top shares

Consider some ordered income $\{y_1,y_2,\dots,y_n\}$, with $y_1\leq y_2\leq\dots\leq y_n$. A classical tool to visualize inequality is Lorenz curve: define the proportion of people $F_{i}=i/n$ (with the convention $F_{0}=0$); then the cumulated wealth $S_{i}=\sum_{j=1}^{i}y_{j}$ and the fraction of cumulated wealth $L_{i}=S_{i}/S_{n}$ (with again ${\displaystyle L_{0}=0}$). Then Lorenz curve is simply the plot $\{F_i,L_i\}$ : it plots the proportion of the total income of the population ($y$ axis) that is cumulatively earned by the bottom $x$\% of the population. And Gini index is the ratio of the area that lies between the line of equality (the first diagonal, $(0,0)-(1,1)$) and the Lorenz curve over the total area under the line of equality. A simple formula would be $${\displaystyle G={\frac {2\sum _{i=1}^{n}iy_{i}}{n\sum _{i=1}^{n}y_{i}}}-{\frac {n+1}{n}}}$$but let us keep in mind simply the fact that it is simply the area below the first diagonal. Note further that the Lorenz curve is increasing, and convex. So actually, for a given Gini index – say $G=60\%$, we can have the two following situations below : on the left, 60% of the poor people get absolutely nothing, and the top 40% shares equally the remaining wealth; on the right, one person gets 60% of the wealth, and everyone else shares equally the remaining wealth.

The two areas are equals (the triangles are the same – up to some symmetrys and rotations) so the two Lorenz curve exhibit the same Gini index. On the left, the 10% the poorest own 0% of the wealth (in green) while the 10% of the richest own 25% of the wealth (in red). On the right the 10% the poorest own 4% of the wealth (in green) while the 10% of the richest own 64% of the wealth (in red). Which can be seen as some sort of paradox : the two cases exhibit the same over inequality, but the one where the poorest get more is also the one where the richest get more.

# More on Random dollars for everyone !

Following my post of yesterday evening, Alex (@AlexSablay) suggested me to look at the Boltzman-Gibbs distribution (e.g. in Yakovenko & Rosser (2009)). There are indeed interesting ideas, and it looks it is more or less what we tried to do in our previous post

Again, I found that article hard to read, but at some point, it looks like they mention that the limiting distribution could be a discrete version that tends to the exponential distribution when the size of the population tends to infinity. Here we have 2000 people, so it should be possible to see it..

If we go for 100,000 rounds, the range of wealth is

so it is still hard to say about the upper bound… For the distribution of the wealth, at the end we obtain the following histogram

and the empirical cumulative distribution function is

Here the red line is the exponential distribution…

So, indeed, it seems that there is a limiting distribution, and it is the exponential one… And the good thing with stable distributions is that they are some sort of fixed point : if we start with that distribution, we should not move (too much) from is. For instance, if we start with an exponential distribution

x = rexp(n,1/init) x = x*init/mean(round(x)) x = round(x)

the range of the wealth remains very stable

as well as the density (again, it is a (symmetric)-kernel based estimate, with a multiplicative bias in 0, and some negative values)

If we plot Lorenz curve, we can see that inequalities do not change here

In that case, it is well known that the Lorenz curve is $u\mapsto u+(1-u)\log(1-u)$ and Gini coefficient is exactly $1/2$.

# Random dollars for everyone !

During the week-end, Philippe Rivière made me discover an interesting problem,

Everyone in a room keeps giving dollars to random others.
You’ll never guess what happens next.f

It was coming from a post, a few years ago on decisionsciencenews.com… This problem was mentioned in recent post since it is related to an article published in the American Scientist in november 2019, Is Inequality Inevitable? (that was translated in French last week, for Pour la Science in a section wrongly entitled Economics since it is only a physicist vision of an (old) economic problem) – see also Brewster Kahle’s post.

(for those really interested in mathematics of inequalities, with a (mathematical) economic perspective, there are countless interesting articles…. see at least Thony Atkinson‘s book or several articles published in Econometrica – references are given in the slides of the course I gave a few years ago on that topic).

I wanted to try, on my own, because I did not understood most of the posts. Because my first thought is that the problem is ill-posed. First of all, what is this “giving dollars”? is it a fixed amount or a random one ? Let us start by assuming that it is fixed. Now, if you know a little bit about gambling and ruin, you guess that it’s very likelely that some one will get banckrupt (at least on a very very long range)… what should we do with that person? Actually, those points were clarified in Jordan’s post

“Imagine a room full of 100 people with 100 dollars each. With every tick of the clock, every person with money gives a dollar to one randomly chosen other person. After some time progresses, how will the money be distributed?”

A well-posed problem states that only people with money can give (everyone can receive) and the amount of money given is fixed.

• A first model (with possible bankruptcy)

First of all, assume that everyone has a fixed amont of money, say 100 (as discussed above). And that each one must give 1 to someone, picked randomly, or more precisely

“every person gives a dollar to one randomly chosen other person”

So, the other people of person $i$ means sampling in $\{1,2,\cdots,n\}\backslash\{i\}$

n = 2000 ns = 20000 init = 100 x = rep(init,n) VX = x VR = c(x[1],x[1]) for(s in 1:ns){ r = function(i) sample((1:n)[-i],size=1) other = Vectorize(r)(1:n) dx = table(other) dx = as.numeric(dx[as.character(1:n)]) dx[is.na(dx)]=0 x = x -rep(1,n)+dx VR=cbind(VR,range(x)) if(s %% 200 ==0) VX=cbind(VX,x) }

Here, I store the range of the wealth of my 2000 people, and every 200 rounds, I also keep tracks of the wealths. The plot of the evolution of the range is the following,

As expected, some people will be ruined… and so far, I did nothing, they keep playing… An easy solution would have been to given them an initial endowment of 1000, and not 100. But that’s only a temporary solution: over 20,000 rounds, there might have no bankruptcy, but over 200,000 there will be ! Before moving to the reflected problem (where only people with money give a dollar), just look at the evolution of the distribution of wealths,

or the evolution of the cumulative distribution

We clearly have more variability as we play. Here, I cannot compute any inequality indices (Lorenz curve is constructed only for positive wealths for instance).

I did not look at analytical results here. The only thing that I know for sure is that about (if there are enough people sharing money) one third (actually $36.78\%$ i.e. $e^{-1}$) will give one dollar, and receive nothing… that’s the law of small numbers (that result was mentioned in Jordan’s post).

• The reflected problem (with no bankruptcy)

Consider now the reflected problem

“Imagine a room full of people with the same amount of money. With every tick of the clock, every person with money gives a dollar to one randomly chosen other person. After some time progresses, how will the money be distributed?”

(I call that reflected because if someone hits the zero-barrier, it can only go up : that person gives nothing, and can possibly receive)

for(s in 1:ns){ r = function(i) sample((1:n)[-i],size=1) other = Vectorize(r)(which(x>0)) dx = table(other) dx = as.numeric(dx[as.character(1:n)]) dx[is.na(dx)] = 0 x = x -(x>0)*1+dx VR = cbind(VR,range(x)) if(s %% 200 ==0) VX = cbind(VX,x) }

Here the range is the following

We are bounded from below (it is not possible to have less than 0) and it seems that extremely reach people are less rich than before. We can look now at the cumulative distribution function (since there is no density here, because of the mass at 0)

(for to get some smooth function, I used a symmetric kernel estimate here, so numerically there are values below 0). Since wealths are positive, we can look at Lorenz curves

It seems that there are more and more inequality, as we play that reallocation game. But here again, I will have to run more simulations (and actually a lot more*) to see if there is a non-degenerated limit with such a game. Here, the distribution of wealth after $n$ rounds is an homogenous Markov chain, taking values in $\mathbb{N}_+$, and using combinatorials, it should be possible to get the transition matrix…

* in did try (during the night) following the advise of Alex (@AlexSablay) advise, and indeed, there is a limiting distribution, see here

• When the contribution is a fixed part (e.g. 1%) of the wealth

An important issue previously was about additivity : “every person with money gives a dollar“. Inequality measures do not like additive operations, they like multiplicate operations (see Serge Christophe Kohlm’s discussion, for instance), or using other words, changes should be relative, not absolute. What about the following question

“Imagine a room full of 100 people with the same amount of money. With every tick of the clock, every person gives a fixed percentage of his money to one randomly chosen other person. After some time progresses, how will the money be distributed?”

The code will be the following: as previously, we match givers and receivers, but here, we have to compute how people give (here it is 1/100 of the money, at each round). At the very first round, we are strictly equivalent to the previous versions : everyone gives 1. The only thing is that, at the second round, those who got nothing at the first one are required to give “only” 99¢.

frac = 1/100 for(s in 1:ns){ r = function(i) sample((1:n)[-i],size=1) other = Vectorize(r)(1:n) df = data.frame(dep = 1:n, arr = other, mont = x*frac) A = aggregate(df$mont,by=list(df$arr),FUN=sum) dx = A$x names(dx) = as.character(A$Group.1) dx = as.numeric(dx[as.character(1:n)]) dx[is.na(dx)] = 0 x = x*(1-frac)+dx VR = cbind(VR,range(x)) if(s %% 200 ==0) VX = cbind(VX,x) }

Here is looks like we have some sort of convergence… at least, no one gets less than 75, and more than 125… The distribution can be visualized below

or via the cumulative distribution function

But to be honest, I don’t know what that distribution is…

To conclude, we can also try something (slightly) different : what if we start with non identical wealths ? Instead of having everyone with wealth 100$, what if it was uniformely distributed between 0$ and 200\$ ?

x = seq(0,2*init,length=n)

It looks like we have a convergence towards the same distribution, with clearly less inequality than when we started… Here is the cumulative distribution (that started with the uniform distribution)

Again, if someone know what that limiting distribution is, I’d be glad to know !

# Pareto models for risk management

Our paper, with Emmanuel Flachaire, “Pareto models for risk management” is now online…

The Pareto model is very popular in risk management, since simple analytical formulas can be derived for financial downside risk measures (Value-at-Risk, Expected Shortfall) or reinsurance premiums and related quantities (Large Claim Index, Return Period). Nevertheless, in practice, distributions are (strictly) Pareto only in the tails, above (possible very) large threshold. Therefore, it could be interesting to take into account second order behavior to provide a better fit. In this article, we present how to go from a strict Pareto model to Pareto-type distributions. We discuss inference, and derive formulas for various measures and indices, and finally provide applications on insurance losses and financial risks.

# Des modeles prédictifs en assurance

Cet post est aussi en ligne sur https://hal.archives-ouvertes.fr/hal-02350006, il a été coécrit avec Laurence Barry et Ewen Gallic.

Les compagnies d’assurance émettent des contrats qui prévoient des paiements d’indemnités en cas de survenance d’évènements aléatoires (accident, maladie, décès, etc.). En contrepartie, l’assuré doit s’acquitter d’une prime, dont le montant est déterminé ex-ante, avant le début de la période de couverture. Cette prime se décompose en deux termes : une prime pure (destinée à couvrir les pertes anticipées) et un chargement (incluant des commissions à des agents, divers frais, mais aussi couvrant contre le risque de variabilité des pertes). La prime pure est souvent calculée par classe de risque, et une classification est alors nécessaire.

# Assurer une population hétérogène, ou l’importance de la classification

Le regroupement des risques selon diverses informations telles l’âge de l’assuré, son état de santé ou encore sa profession constitue ce que l’on appelle la classification des risques. Cette pratique de segmentation se justifie (à des fins d’admissibilité mais aussi de tarification) par la supposition que les risques sont placés dans des groupes relativement homogènes, au sein desquels les probabilités de survenance sont similaires. Pour Schauer (2006), cette « généralisation », qui vise à voir l’individu sous le prisme de sa classe de risque, de généraliser son comportement à partir de quelques variables explicatives, est probablement la raison d’être de l’actuaire : « to be an actuary is to be a specialist  in  generalization,  and  actuaries  engage  in  a form of decision-making that is sometimes called actuarial ». Statistiquement on cherche une méthode de classification aussi « discriminatoire » que possible[1], en gardant en mémoire que la discrimination est interdite, ce qui rend l’exercice périlleux, et souvent critiqué (nous y reviendrons plus loin).

Les assureurs évoquent souvent deux arguments pour justifier une segmentation. Le premier est qu’elle serait rendue économiquement nécessaire par la concurrence ; ne pas classifier conduit à une anti-sélection, les risques importants restant seuls chez les assureurs qui ne segmentent pas. Dans une telle situation, l’équilibre de marché ne serait pas possible puisque les risques faibles seraient chez un concurrent ayant segmenté. Si le facteur de risque était observable, tant par les assurés que les assureurs, il y aurait un phénomène d’auto-sélection, les assurés à risque faible ayant les polices les moins chères. Cette situation constitue un équilibre séparant de Nash. Mais si le facteur de risque n’est pas observable, un équilibre sous-optimal peut être atteint, résultant d’une externalité négative de cette information non-accessible, à la manière de Wilson (1977), tel que décrit dans Cummins et al. (1982) dans le cas des contrats d’assurance-vie. Cela dit, Kleindorfer & Kunreuther (1980) montrent qu’accéder à davantage d’information ne conduit pas nécessairement à une amélioration du bien-être des consommateurs. De plus si la classification n’est pas autorisée, l’équilibre est maintenu, les risques faibles subventionnant les risques élevés.

Le second argument avancé pour justifier une segmentation est que cette dernière (et le fait, par conséquent, d’ajuster les primes au risque) serait juste et équitable. Mais cette vision de l’équité n’a pas toujours été de mise et semble portée par les développements techniques. Ainsi la classification est devenue de plus en plus fine, multipliant les classes de risque et conduisant à des tarifs « personnalisés ». En plus des avancées statistiques, des facteurs économiques pourraient justifier cette sophistication : la concurrence de plus en plus forte sur certaines branches.

# Incertitude en assurance

Il y a plusieurs manières de caractériser l’incertitude en assurance. Comme bien souvent quand on fait des prévisions, il convient de distinguer l’incertitude associée à l’estimation des probabilités et l’incertitude réelle sur le résultat (aléa de l’évènement). Pour la seconde notion, Hacking (1975) parle de probabilité structurelle, et c’est celle qui est souvent utilisée pour introduire les concepts de probabilité, par exemple avec des dés ou des jeux de cartes : les probabilités sont connues, seule l’issue du jeu est incertaine. Par exemple je sais que la probabilité d’avoir 6 en lançant un dé est 1/6 (compte tenu de la géométrie du cube).

D’un point de vue statistique, la probabilité se mesure quand on peut observer une fréquence, c’est-à-dire une répétition de risques semblables. Les statisticiens ont ainsi défini une notion de probabilité empirique, basée sur la répétition[2]. Si, en lançant mille dés j’obtiens 173 fois la face 6, la probabilité empirique d’avoir 6 est de 17,3%. La loi des grands nombres nous assure que cette fréquence va tendre vers la vraie valeur en répétant l’expérience, et le théorème central limite permet d’en contrôler les fluctuations. C’est la première incertitude dont nous parlions au début de cette section, que nous appellerions l’erreur d’estimation.

On peut enfin mentionner deux notions supplémentaires ; tout d’abord, les probabilités conditionnelles. Cette idée est introduite en assurance par de Moivre, ou de Witt, lorsqu’ils notaient que pour estimer une probabilité de décès, il fallait considérer des personnes de même âge. C’est cette idée que l’on retrouve quand on considère une classification : on veut des risques homogènes, similaires, sans être pour autant identiques. La probabilité que l’on obtient est alors conditionnelle à ce facteur commun qui caractérise la classe observée. Dans notre exemple des dés, cela revient à dire qu’il ne faut pas lancer mille dés, mais mille fois le même dé – ou à défaut des dés semblables.

Enfin, les probabilités subjectives ont été formalisées par Bruno de Finetti et Leonard Savage (ainsi que plus philosophiquement par Frank P. Ramsey) pour comprendre et modéliser la prise de décision. Elles sont relativement populaires en économie de l’incertain, mais difficile à mettre en œuvre dans un contexte de valorisation de contrats d’assurance automobile ou habitation. Il s’agit d’un jugement, qui ne peut être confronté à la réalité, mais envisageable pour l’assurance de risques encore mal connus (McGrayne (2012) évoque ainsi les premiers contrats d’assurance aviation). Une approche bayésienne consiste alors à combiner cette probabilité subjective avec la probabilité comme fréquence observée d’un phénomène : partant d’une croyance a priori, on affine l’estimation par une mise à jour progressive en répétant les expériences. Classiquement, la probabilité d’avoir la face 6 sera une moyenne entre notre croyance (1 chance sur 6) et une probabilité dite historique, obtenue en faisant quelques lancés (3 sur 20 lancers, par exemple). Les poids attribués aux deux dépendant du nombre d’expériences effectuées : on donnera plus de crédit à l’expérience si on fait mille lancés que si on en fait soixante.

# Incertitude sur le résultat, ou aléa fondamental

Les probabilités prédictives, utilisées pour calculer la prime d’un contrat d’assurance, sont la première étape d’un problème de classification. Un outil classique pour juger de la pertinence d’un classifieur est la courbe ROC, décrite dans Kuhn (2018)) : on compare la probabilité individuelle (a priori, telle que résultant du modèle de classification) à un seuil, compris entre 0 et 1; si la probabilité est inférieure au seuil, l’estimation est que la personne survit, sinon qu’elle décède.

On compare ensuite cette estimation aux réalisations (ex-post) de survie et de décès. Pour chaque seuil, on peut considérer la matrice classique dite matrice de confusion de théorie de la décision : elle consiste à répartir les observations suivant le résultat observé (en colonne) et l’estimation résultant du modèle en ligne (en fonction de la probabilité estimée pour l’individu et le seuil que l’on s’est fixé). On peut ainsi partager la population entre les classements corrects, et les erreurs (dont les « faux positifs » si la personne a survécu malgré une probabilité estimée de décès supérieure au seuil, et les « faux négatifs » si la personne décède malgré une probabilité estimée inférieure au seuil).

Figure 1: Courbe ROC et classification pour un seuil de probabilité valant 1.5%.

La courbe ROC est obtenue en faisant varier le seuil. Chaque seuil correspond à un point de la courbe, rapportant graphiquement les taux de faux positifs (en abscisse) et de vrais positifs (en ordonnée), comme sur la Figure 1.

Considérons un groupe de 1000 assurés, où 20 personnes sont décédées l’an passé. Supposons un modèle dans lequel on admet que la population est parfaitement homogène, la probabilité estimée de décès est de 2% pour tout le monde. Dans ce cas pour tout seuil supérieur à 2%, on estimera que la totalité de la population survit : on aura un taux de faux positifs de 0% et un taux de vrais positifs de 0%, d’où un point (0,0) sur le graphe. A l’inverse pour tout seuil inférieur à 2%, on estimera que la totalité de la population décède : on aura un taux de faux positifs de 100% et un taux de vrais positifs de 100%, d’où un point (1,1) sur le graphe. La courbe de ROC de ce modèle uniforme à 2% est donc la diagonale du carré sur la figure 1.

Mais on peut aussi imaginer qu’il existe un peu d’hétérogénéité avec, par exemple, une probabilité de décès de 1% pour une moitié de la population et de 3% pour l’autre moitié, ou encore que le modèle produit des probabilités comprises entre 1% et 3% de façon non dichotomique. Les données simulées pour construire la courbe noire sur la Figure 1 suppose que la population a des probabilités de décès variables, comprises entre 1% et 3%, obtenues par une régression logistique. Comme le montre le tableau de droite, on commet des erreurs, et comme le montre celle de gauche, la nature de celle-ci varie en fonction du seuil choisi, qui modifie les taux de faux positifs et de faux négatifs.

Le cas extrême serait celui où le modèle aurait correctement attribué une probabilité de 100% aux 20 personnes qui sont effectivement décédées. C’est la courbe rouge sur la Figure 1. Ce partage est possible ex-post, une fois réalisation de l’aléa : a posteriori, il y a une certitude de décès pour ceux qui sont effectivement morts. Mais cela n’a cependant pas grande réalité dans l’assurance, à moins d’imaginer que l’actuaire serait un oracle, qui saurait avec certitude qui va mourir, et qui va survivre. La réalité est plutôt celle de la situation intermédiaire entre la courbe rouge et la diagonale, avant d’arriver dans la région hachurée, où le taux d’erreur est faible, mais pas nul : on ne peut pas prédire, avec certitude, qui va décéder. L’assurance n’est possible que si cette borne supérieure n’est pas trop élevée.

# Incertitude statistique, données et modèles

Une question fondamentale pour la survie de l’assurance est de savoir où se situe cette borne supérieure : jusqu’où peut-on aller, entre les deux cas extrêmes (population homogène avec une probabilité de 2% pour tous, et une population très discriminée, avec 2% de la population ayant 100% de chances de mourir, et l’autre 0%) ? Et de quoi cette borne dépend-elle ? En particulier, des modèles plus complexes, tels que les réseaux de neurones très profonds permettent-ils vraiment d’améliorer la prévision ? Et l’enrichissement de données, tel qu’on l’observe grâce aux objets connectés et la fusion avec toutes sortes d’informations externes, va-t-il déplacer la borne supérieure vers le haut ?

Si l’apprentissage profond – voir Goodfellow et al. (2018) – permet d’avoir des classifieurs d’images avec un taux d’erreur proche de 0%, il est difficile d’imaginer qu’il sera possible de prévoir, presque un an à l’avance (à la signature du contrat), qui décèdera dans l’année, qui aura la grippe, qui aura un dégât des eaux, etc. Les modèles plus complexes permettent d’améliorer les prévisions, en tenant compte de non-linéarités, d’effets croisés entre les variables tarifaires, mais pas au point de faire disparaître l’aléa. Et tant que l’assurance est envisagée ex-ante (la prime est fixée au début de la période de couverture), il est difficile d’imaginer que rajouter de l’information fera aussi disparaître l’aléa. C’est d’ailleurs le cas pour les tests génétiques qui n’expliquent qu’une (petite) partie du risque de cancer, par exemple. Et rajouter des données revient souvent à rajouter du bruit, ce qui rend le travail d’analyse plus complexe. Cependant, force est de constater que des modèles plus complexes et des données plus riches tendent effectivement à « améliorer » la prévision, en remontant la courbe ROC vers le haut. Mais se pose-t-on les bonnes questions ? Que signifie vraiment une borne très éloignée du cas homogène, sur la diagonale ?

# Homogénéité, équité et causalité

Comme nous l’avons vu, la tarification en assurance repose sur une répartition des risques (des contrats) en catégories, au sein desquelles la distribution des pertes peut être estimée, afin de fixer un niveau de primes. La répartition se fait à partir des caractéristiques de l’assuré, et du bien assuré. En retraçant l’histoire de l’assurance, Ewald (1986) montre que les mécanismes de prévoyance se sont mis en place en déplaçant la charge des accidents du travail sur la société : on abandonne l’idée d’une responsabilité individuelle de l’accident en faveur de la solidarité. L’assurance distingue « entre le dommage que subit tel ou tel individu — c’est affaire de chance ou de malchance — et la perte liée au dommage dont l’attribution est, quant à elle, toujours collective et sociale ». Ce principe de solidarité sociale, de mutualisation des risques, fait que le risque (en assurance) est toujours pensé collectivement.

Aujourd’hui, les tarifs sont considérés comme « justes », ou « actuariellement équitables » si chaque prime correspond à la perte attendue (pour ne pas dire « espérée », au sens mathématique) pour chaque assuré. Dans cette perception de l’équité, une hypothèse essentielle est que les classes soient « homogènes ». En effet, dans l’hypothèse inverse, les personnes les moins risquées subventionnent les personnes les plus risquées, ce qui est perçu comme socialement injuste.

On peut décrire cette version de l’équité actuarielle à l’aide de la formule de décomposition de la variance. La variance globale se décompose en effet en deux termes, la variance inter-classes et la variance intra-classes : l’ « équité actuarielle » vise à ce que  les classes de risque soient relativement distinctes les unes des autres, donc une variance inter-classes forte, accompagnée d’une homogénéité des classes, donc une variance intra-classes faible. D’un point de vue statistique, chercher à augmenter l’une est équivalente à faire diminuer l’autre. Cette mécanique n’est pas toujours claire pour des observateurs non avertis ; ainsi dans l’affaire Manhart, un des cas les plus documentés sur la discrimination par le genre en assurance, le juge Stevens affirme : « we focus on fairness to individuals rather than on fairness to classes […] even a true generalization about a class is an insufficient reason for disqualifying an individual to whom the generalization does not apply» (cité dans Anzalone (2016)). Autrement dit, pour la justice, un critère statistique de type « true generalization » ne peut être opposable à un individu.

Une autre critique importante, que l’on retrouve dans la « gender directive », est le lien entre discrimination et causalité. En effet, statistiquement, les actuaires vont chercher des facteurs de classification fortement corrélés avec la sinistralité. Mais il est possible que ces facteurs ne soient qu’un proxy de la vraie variable causale, restée elle inobservée, conduisant à une mauvaise estimation du risque pour certains. Comme le notent Antonio et Charpentier (2017), le genre a ainsi été longtemps utilisé en assurance automobile car très corrélé avec des variables associées au style de conduite et à d’autres variables historiquement non observables (mais qui le sont aujourd’hui grâce aux objets connectés, comme le kilométrage, les heures de conduites, les types de routes utilisés, etc).

Ce lien avec les mécanismes causaux est d’ailleurs relativement profond, et Hacking (1975) y voit une connexion avec la « révolution probabiliste » : on peut assez facilement mettre en évidence des corrélations, mais les causes, si elles existent, nous restent plus opaques. Laplace au début du 19e siècle déclare ainsi que « la probabilité est relative en partie à nos connaissances, en partie à notre ignorance », liant les probabilités à la fois à une vision newtonienne déterministe du monde et à notre incapacité à le connaitre parfaitement. Cette dernière composante fait que l’on ne peut pas annoncer la date exacte du décès d’un individu, mais statistiquement, dans un groupe homogène, on peut prédire le nombre de décès au cours d’une année. Et pour revenir à la relation causale, le tabagisme par exemple ne cause pas forcément une mort prématurée mais fumer sera vu comme dangereux car il augmente la probabilité de décès pendant une période donnée. Ainsi nous montre Hacking (1975), la causalité est pensée aujourd’hui dans un contexte probabiliste, et non plus déterministe.

Abraham, K.S. (1985) Efficiency and fairness in insurance risk classification, Virginia Law Review 71: 403-451.

Anderson, A.W. (1978). A Critique of the Manhart Brief. The Actuary, 12:5.

Antonio, K. & Charpentier, A. (2017). La tarification par genre en assurance, corrélation ou causalité ? Risques, 109.

Anzalone, C.A. (2016). U.S. Supreme Court Cases on Gender and Sexual Equality. Routledge.

Bailey, H., Hutchison, T. & Narber, G. (1975) The regulatory challenge to life insurance classification, Drake Law Review Insurance Law Annual 4: 779-827

Barry L. (2019). Justice ou justesse? L’équité de l’assurance. Working paper #15, chaire PARI.

Charpentier, A. & Denuit, M. (2004). Mathématiques de l’Assurance Non-Vie : Principes Généraux de Théorie du Risque. Economica.

Cramer, H. (1946). Mathematical Methods of Statistics. Princeton University Press.

Cummins, J.D., Smith, B.D., Vance, R.N. & VanDerhai, J.L. (1982). Risk Classification in Life Insurance. Kluwer-Nijhoff Publishing.

Ewald F. (1986). L’État providence. Grasset.

Fisher, R. A. (1936).  The Use of Multiple Measurements in Taxonomic Problems. Annals of Eugenics. 7 (2): 179–188.

Frézal, S. & Barry, L. (2019). Fairness in Uncertainty: Some Limits and Misinterpretations of Actuarial Fairness, Journal of Business Ethics.

Goodfellow, U., Bengio, Y. & Courville, A. (2018) L’apprentissage profond. Massot Edition.

Hacking, I. (1975) The Emergence of Probability. Cambridge University Press

Hoy, M. (1982) Categorizing risks in the insurance industry, The Quarterly Journal of Economics 97: 321-336.

Kleindorfer, P. & Kunreuther, H. (1980) Misinformation and Equilibrium in Insurance Markets, in Economic Analysis of Regulated Markets, Jörg Finsinger Editor, Springer Verlag, 67-90

Kuhn, M. & Johnson, K. (2018). Applied Predictive Modeling. Springer Verlag.

Martin,G.D. (1976) Gender Discrimination in Pension Plans, Journal of Risk and Insurance 43.

McGrayne, S.B. (2012) The Theory That Would Not Die: How Bayes’ Rule Cracked the Enigma Code, Hunted Down Russian Submarines, and Emerged Triumphant from Two Centuries of Controversy. Yale University Press.

Ramsey P.F. (1926). Truth and Probability.

Schauer, F. (2006) Profiles, Probabilities, and Stereotypes. Harvard University Press.

Von Mises, R. (1957). Probability, Statistics and Truth. Dover publications.

Wilson, C. (1977). A model of insurance markets with incomplete information. Journal of Economic Theory, 16:2, 167-207.

[1] Au sens statistique du mot, dans le sens introduit par Fisher (1936).

[2] Dans cette approche fréquentiste, et notamment pour Ronald Fisher et Richard von Mises, la probabilité d’un évènement unique (dit « one shot ») n’a pas de sens.

# Big Data, GAFA et Assurance

Les sociétés technologiques et le monde de l’assurance auraient tout pour être opposé. Agilité, rapidité, obsession du futur chez les uns, conservatisme, réflexivité, fascination pour les données passées chez les autres. Et pourtant les deux s’observent, et commencent à nouer des partenariats, comprenant que la donnée est leur cœur de métier.

# Foundations of Machine Learning, part 5

This post is the nineth (and probably last) one of our series on the history and foundations of econometric and machine learning models. The first fours were on econometrics techniques. Part 8 is online here.

## Optimization and algorithmic aspects

In econometrics, (numerical) optimization became omnipresent as soon as we left the Gaussian model. We briefly mentioned it in the section on the exponential family, and the use of the Fisher score (gradient descent) to solve the first order condition $\mathbf{X}^T W(\beta)^{-1})[y-\widehat{y}]=\mathbf{0}$. In learning, optimization is the central tool. And it is necessary to have effective optimization algorithms, to solve problems (described previously) of the form: $$\widehat{\beta}\in\underset{\beta\in\mathbb{R}^p}{\text{argmin}}\left\lbrace\sum_{i=1}^n \ell(y_i,\beta_0+\mathbf{x}^T\beta)+\lambda\Vert\boldsymbol{\beta}\Vert\right\rbrace$$In some cases, instead of global optimization, it is sufficient to consider optimization by coordinates (widely studied in Daubechies et al. (2004)). If $f:\mathbb{R}^d\rightarrow\mathbf{R}$ is convex and differentiable, if $\mathbf{x}$ satisfies $f(\mathbf{x}+h\boldsymbol{e}_i)\geq f(\mathbf{x})$ for any $h>0$ and $i\in\{1,\cdots, d\}$then $f(\mathbf{x})=\min\{f\}$, where $\mathbf{e}=(\mathbf{e}_i)$ is the canonical basis of $\mathbb{R}^d$. However, this property is not true in the non-differentiable case. But if we assume that the non-differentiable part is separable (additively), it becomes true again. More specifically, if$$f(\mathbf{x})=g(\mathbf{x})+\sum_{i=1}^d h_i(x_i)$$with$$\left\lbrace\begin{array}{l}g: \mathbb{R}^d\rightarrow\mathbb{R}\text{ convex-differentiable}\\h_i: \mathbb{R}\rightarrow\mathbb{R}\text{ convex}\end{array}\right.$$This was the case for Lasso regression, $\beta)\mapsto\| \mathbf{y}-\beta_0-\mathbf{X}\beta\|_{\ell_2 }+\lambda\|\beta\|_{\ell_1}$, as shown by Tsen (2001). Getting back to our initial notations, we can use a coordinate descent algorithm: from an initial value $\mathbf{x}^{(0)}$, we consider (by iterating)$$x_j^{(k)}\in\text{argmin}\big\lbrace f(x_1^{(k)},\cdots,x_{k-1}^{(k)},x_k,x_{k+1}^{(k-1)},\cdots,x_n^{(k-1)})\big\rbrace$$ for $j=1,2,\cdots,n$These algorithmic problems and numerical issues may seem secondary to econometricians. However, they are essential in automatic learning: a technique is interesting if there is a stable and fast algorithm, which allows to obtain a solution. These optimization techniques can be transposed: for example, this coordinate descent technique can be used in the case of SVM methods (known as “vector support” methods) when the space is not linearly separable, and the classification error must be penalized (we will come back to this technique in the next section).

## In-sample, out-of-sample and cross-validation

These techniques seem intellectually interesting, but we have not yet discussed the choice of the penalty parameter $\lambda$. But this problem is actually more general, because comparing two parameters $\widehat{\beta}_{\lambda_1}$ and $\widehat{\beta}_{\lambda_2}$ is actually comparing two models. In particular, if we use a Lasso method, with different thresholds $\lambda$, we compare models that do not have the same dimension. Previously, we have addressed the problem of model comparison from an econometric perspective (by penalizing overly complex models). In the learning literature, judging the quality of a model on the data used to construct it does not make it possible to know how the model will behave on new data. This is the so-called “generalization” problem. The traditional approach then consists in separating the sample (size $n$) into two parts: a part that will be used to train the model (the training database, in-sample, size $m$) and a part that will be used to test the model (the testing database, out-of-sample, size $n-m$). The latter then makes it possible to measure a real predictive risk. Suppose that the data are generated by a linear model $y_i=\mathbf{x}_i^T \beta_0+\varepsilon_i$ where $\varepsilon_i$ are independent and centred law achievements. The empirical quadratic risk in-sample is here$$\frac{1}{m}\sum_{i=1}^m\mathbb{E}\big([\mathbf{x}_i^T \widehat{\beta}-\mathbf{x}_i^T \beta_0]^2\big)=\mathbb{E}\big([\mathbf{x}_i^T \widehat{\beta}-\mathbf{x}_i^T \beta_0]^2\big),$$for any observation $i$. Assuming the residuals $\varepsilon$ Gaussian, then we can show that this risk is worth $\sigma^2 \text{trace} (\Pi_X)/m$ is $\sigma^2 p/m$. On the other hand, the empirical out-of-sample quadratic risk is here $$\mathbb{E}\big([\mathbf{x}^T \widehat{\beta}-\mathbf{x}^T \beta_0]^2\big)$$where $\mathbf{x}$ is a new observation, independent of the others. It can be noted that $$\mathbb{E}\big([\mathbf{x}^T \widehat{\beta}-\mathbf{x}^T \beta_0]^2\big\vert \mathbf{x}\big)=\text{Var}\big(\mathbf{x}^T \widehat{\beta}\big\vert \mathbf{x}\big)=\sigma^2\mathbf{x}^T(\mathbf{x}^T\mathbf{x})^{-1}\mathbf{x},$$and by integrating with respect to $\mathbf{x}$, $$\mathbb{E}\big([\mathbf{x}^T \widehat{\beta}-\mathbf{x}^T\beta_0]^2\big)=\sigma^2\text{trace}\big(\mathbb{E}[\mathbf{x}\mathbf{x}^T]\mathbb{E}\big[(\mathbf{x}^T\mathbf{x})^{-1}\big]\big).$$The expression is then different from that obtained in-sample, and using the Groves & Rothenberg (1969) increase, we can show that $$\mathbb{E}\big([\mathbf{x}^T \widehat{\beta}-\mathbf{x}^T \beta_0]^2\big) \geq \sigma^2\frac{p}{m}$$which is pretty intuitive, when we start thinking about it. Except in some simple cases, there is no simple (explicit) formula. Note, however, that if $\mathbf{X}\sim\mathcal{N}(0,\sigma^2 \mathbb{I})$, then $\mathbf{x}^T \mathbf{x}$ follows a Wishart law, and it can be shown that $$\mathbb{E}\big([\mathbf{x}^T \widehat{\beta}-\mathbf{x}^T \beta_0]^2\big)=\sigma^2\frac{p}{m-p-1}.$$If we now look at the empirical version: if $\widehat{\beta}$ is estimated on the first $m$ observations,$$\widehat{\mathcal{R}}^{~\text{ IS}}=\sum_{i=1}^m [y_i-\boldsymbol{x}_i^T\widehat{\boldsymbol{\beta}}]^2\text{ and }\widehat{\mathcal{R}}^{\text{ OS}}=\sum_{i=m+1}^{n} [y_i-\boldsymbol{x}_i^T\widehat{\boldsymbol{\beta}}]^2$$and as Leeb (2008) noted, $\widehat{\mathcal{R}}^{\text{IS}}-\widehat{\mathcal{R}}^{\text{OS}}\approx 2\cdot\nu$ where $\nu$ represents the number of degrees of freedom, which is not unlike the penalty used in the Akaike test.

Figure 4 shows the respective evolution of $\widehat{\mathcal{R}}^{\text{IS}}$ and $\widehat{\mathcal{R}}^{\text{OS}}$ according to the complexity of the model (number of degrees in a polynomial regression, number of nodes in splines, etc). The more complex the model, the more $\widehat{\mathcal{R}}^{\text{IS}}$ will decrease (this is the red curve, below). But that’s not what we’re interested in here: we want a model that predicts well on new data (i. e. out-of-sample). As Figure 4 shows, if the model is too simple, it does not predict well (as it does with in-sample data). But what we can see is that if the model is too complex, we are in a situation of “overlearning”: the model will start to model the noise. Of course, this figure should remind us of the one we’ve seen in our second post of that series

Figure 4 : Generalization, under- and over-fitting

Instead of splitting the database in two, with some of the data that will be used to calibrate the model and some to study its performance, it is also possible to use cross-validation. To present the general idea, we can go back to the “jackknife”, introduced by Quenouille (1949) (and formalized by Quenouille (1956) and Tukey (1958)) relatively used in statistics to reduce bias. Indeed, if we assume that $\{y_1,\cdots,y_n\}$ is a sample drawn according to a law $F_\theta$, and that we have an estimator $T_n (\mathbf{y})=T_n (y_1,\cdots,y_n)$, but that this estimator is biased, with $\mathbf{E}[T_n (\mathbf{Y})]=\theta+O(n^{-1})$, it is possible to reduce the bias by considering $$\widetilde{T}_n(\mathbf{y})=\frac{1}{n}\sum_{i=1}^n T_{n-1}(\mathbf{y}_{(i)})\text{ where }\mathbf{y}_{(i)}=(y_1,\cdots,y_{i-1},y_{i+1},\cdots,y_n)$$It can then be shown that $\mathbb{E}[\tilde{T}_n(Y)]=\theta+O(n^{-2})$The idea of cross-validation is based on the idea of building an estimator by removing an observation. Since we want to build a predictive model, we will compare the forecast obtained with the estimated model, and the missing observation$$\widehat{\mathcal{R}}^{\text{ CV}}=\frac{1}{n}\sum_{i=1}^n \ell(y_i,\widehat{m}_{(i)}(\mathbf{x}_i))$$We will speak here of the “leave-one-out” (loocv) method.

This technique reminds us of the traditional method used to find the optimal parameter in exponential smoothing methods for time series. In simple smoothing, we will construct a forecast from a time series as ${}_t\widehat{y}_{t+1} =\alpha\cdot{}_{t-1}\widehat{y}_t +(1-\alpha)\cdot y_t$, where $\alpha\in[0,1]$, and we will consider as “optimal” $$\alpha^\star = \underset{\alpha\in[0,1]}{\text{argmin}}\left\lbrace \sum_{t=2}^T \ell({}_{t-1}\widehat{y}_{t},y_{t}) \right\rbrace$$as described by Hyndman et al (2009).

The main problem with the leave-one-out method is that it requires calibration of n models, which can be problematic in large dimensions. An alternative method is cross validation by $k$-blocks (called “$k$-fold cross validation”) which consists in using a partition of $\{1,\cdots,n\}$ in $k$ groups (or blocks) of the same size, $\mathcal{I}_1,\cdots,\mathcal{I}_k$, and let us note $\mathcal{I}_{\bar j}=\{1,\cdots,n\}\setminus \mathcal{I}_j$. By noting $\widehat{m}_{(j)}$ built on the sample $\mathcal{I}_{\bar j}$, we then set:$$\widehat{\mathcal{R}}^{k-\text{ CV}}=\frac{1}{k}\sum_{j=1}^k \mathcal{R}_j\text{ where }\mathcal{R}_j=\frac{k}{n}\sum_{i\in\mathcal{I}_{{j}}} \ell(y_i,\widehat{m}_{(j)}(\mathbf{x}_i))$$Standard cross-validation, where only one observation is removed each time (loocv), is a special case, with $k=n$. Using $k=5$ or $10$ has a double advantage over $k=n$: (1) the number of estimates to be made is much smaller, 5 or 10 rather than $n$; (2) the samples used for estimation are less similar and therefore less correlated to each other, which tends to avoid excess variance, as recalled by James et al. (2013).

Another alternative is to use boosted samples. Let $\mathcal{I}_b$ be a sample of size $n$ obtained by drawing with replacement in $\{1,\cdots,n\}$ to know which observations $(y_i,\mathbf{x}_i)$ will be kept in the learning population (at each draw). Note $\mathcal{I}_{\bar b}=\{1,\cdots,n\}\setminus\mathcal{I}_b$. By noting $\widehat{m}_{(b)}$ built on sample $\mathcal{I}_b$, we then set :$$\widehat{\mathcal{R}}^{\text{ B}}=\frac{1}{B}\sum_{b=1}^B \mathcal{R}_b\text{ where }\mathcal{R}_b=\frac{n_{\overline{b}}}{n}\sum_{i\in\mathcal{I}_{\overline{b}}} \ell(y_i,\widehat{m}_{(b)}(\mathbf{x}_i))$$where $n_{\bar b}$ is the number of observations that have not been kept in $\mathcal{I}_b$. It should be noted that with this technique, on average $e^{-1}\sim36.7\%$ of the observations do not appear in the boosted sample, and we find an order of magnitude of the proportions used when creating a calibration sample, and a test sample. In fact, as Stone (1977) had shown, the minimization of AIC is to be compared to the cross-validation criterion, and Shao (1997) showed that the minimization of BIC corresponds to $k$-fold cross-validation, with $k=n/\log n$.

All those techniques here are mentioned in the “machine learning” section since they rely on automatic, computational techniques, and no probabilistic foundations are necessary. In many cases we did use the notation $m^\star$ (at least in the first posts on “machine learning” techniques) to highlight the fact that we want some sort of “optimal” model – and to make a distinction with estimators $\widehat{m}$ considered earlier, when we had some probabilistic framework. But of course, it is possible (and necessary) to build bridges between those two cultures…

References are online here. As explained in the introduction, it is some sort of online version of an introduction to our joint paper with Emmanuel Flachaire and Antoine Ly, Econometrics and Machine Learning (initially writen in French), that will actually appear soon in the journal Economics and Statistics (in English and in French).

# Foundations of Machine Learning, part 4

This post is the eighth one of our series on the history and foundations of econometric and machine learning models. The first fours were on econometrics techniques. Part 7 is online here.

## Penalization and variables selection

One important concept in econometrics is Ockham’s razor – also known as the law of parsimony (lex parsimoniae) – which can be related to abductive reasoning.

Akaike’s criterion was based on a penalty of likelihood taking into account the complexity of the model (the number of explanatory variables retained). If in econometrics, it is customary to maximize the likelihood (to build an asymptotically unbiased estimator), and to judge the quality of the ex-post model by penalizing the likelihood, the strategy here will be to penalize ex-ante in the objective function, even if it means building a biased estimator. Typically, we will build: $$(\widehat{\beta}_{0,\lambda},\widehat{\beta}_{\lambda})=\text{argmin}\left\lbrace\sum_{i=1}^n \ell(y_i,\beta_0+\mathbf{x}^T\beta)+\lambda \text{ penalization}( \boldsymbol{\beta})\right\rbrace, ~~~(11)$$where the penalty function will often be a norm $\|\cdot\|$ chosen a priori, and a penalty parameter $\lambda$ (we find in a way the distinction between AIC and BIC if the penalty function is the complexity of the model – the number of explanatory variables retained). In the case of the $\ell_2$ norm, we find the ridge estimator, and for the $\ell_1$ norm, we find the lasso estimator (“Least Absolute Shrinkage and Selection Operator”). The penalty previously used involved the number of degrees of freedom of the model, so it may seem surprising to use $\|\beta\|_{\ell_2}$ as in the ridge regression. However, we can envisage a Bayesian vision of this penalty. It should be recalled that in a Bayesian model : $$\underbrace{\mathbb{P}[\boldsymbol{\theta}\vert\boldsymbol{y}]}_{\text{posterior}} \propto \underbrace{\mathbb{P}[\boldsymbol{y}\vert\boldsymbol{\theta}]}_{\text{likelihood}} \cdot \underbrace{\mathbb{P}[\boldsymbol{\theta}]}_{\text{prior}}$$or$$\log\mathbb{P}[\boldsymbol{\theta}\vert\boldsymbol{y}]= \underbrace{\log \mathbb{P}[\boldsymbol{y}\vert\boldsymbol{\theta}]}_{\text{log likelihood}} + \underbrace{\log\mathbb{P}[\boldsymbol{\theta}]}_{\text{{penalty}}}$$In a Gaussian linear model, if we assume that the a priori law of $\theta$ follows a centred Gaussian distribution, we find a penalty based on a quadratic form of the components of $\theta$.

Before going back in detail to these two estimators, obtained using the $\ell_1$ or $\ell_2$ norm, let us return for a moment to a very similar problem: the best choice of explanatory variables. Classically (and this will be even more true in large dimension), we can have a large number of explanatory variables, $p$, but many are just noise, in the sense that $\beta_j=0$ for a large number of $j$. Let $s$ be the number of (really) relevant covariates, $s=\#S$, with $$S=\{j=1,\cdots,p:\beta_j\neq 0\}$$. If we note $\mathbf{X}_S$ the matrix composed of the relevant variables (in columns), then we assume that the real model is of the form $y=\mathbf{x}_S^T \beta_S+\varepsilon$. Intuitively, an interesting estimator would then be $\widehat{\beta}_S=[\mathbf{X}_S^T \mathbf{X}_S ]^{-1} \mathbf{X}_S^T \mathbf{y}$, but this estimator is only theoretical because the set $S$ is unknown, here. This estimator can actually be seen as the oracle estimator mentioned above. One may then be tempted to solve $$(\widehat{\beta}_{0,s},\widehat{\beta}_{s})=\underset{\beta_S\in\mathbb{R}^s}{\text{argmin}}\left\lbrace\sum_{i=1}^n \ell(y_i,\beta_0+\mathbf{x}^T\beta_S)\right\rbrace,\text{ s.t. } \# {S}=s$$This problem was introduced by Foster & George (1994) using the $\ell_0$ notation. More precisely, let us define here the following three norms, where $\mathbf{a}\in\mathbb{R}^d$, $$\Vert\boldsymbol{a} \Vert_{\ell_0}=\sum_{i=1}^d \mathbf{1}(a_i\neq 0), ~~ \Vert\mathbf{a} \Vert_{\ell_1}=\sum_{i=1}^d |a_i|~~\text{ and }~~\Vert\mathbf{a} \Vert_{\ell_2}=\left(\sum_{i=1}^d a_i^2\right)^{1/2}$$

Table 1: Constrained optimization and regularization.

Let us consider the optimization problems in Table 1. If we consider the classical problem where the quadratic norm is used for $\ell$, the two problems of the equation $(\ell1)$ of Table 1 are equivalent, in the sense that, for any solution $(\beta^\star,s)$ to the left problem, there is $\lambda^\star$ such that $(\beta^\star,\lambda^\star)$ is the solution of the right problem; and vice versa. The result is also true for problems$(\ell2)$. These are indeed convex problems. On the other hand, the two problems $(\ell0)$ are not equivalent: if for $(\beta^\star,\lambda^\star)$ solution of the right problem, there is $s^\star$ such that $\beta^\star$ is solution of the left problem, the reverse is not true. More generally, if you want to use an $\ell_p$ norm, sparsity is obtained if $p\leq 1$ whereas you need $p\geq1$ to have the convexity of the optimization program.

One may be tempted to resolve the penalized program $(\ell0)$ directly, as suggested by Foster & George (1994). Numerically, it is a complex combinatorial problem in large dimension (Natarajan (1995) notes that it is a NP-difficult problem), but it is possible to show that if $\lambda\sim\sigma^2 \log(p)$, then $$\mathbb{E}\big([\mathbf{x}^T \widehat{\beta}-\mathbf{x}^T \beta_0]^2\big) \leq \underbrace{\mathbb{E}\big(\mathbf{x}_{ {S}}^T\widehat{\beta}_{{S}}-\mathbf{x}^T \beta_0]^2\big)}_{=\sigma^2 \#{S}}\cdot \big(4\log p+2+o(1)\big)$$Observe that in this case $$\widehat{\beta}_{\lambda,j}^{\text{sub}} = \left\lbrace\begin{array}{l}0 \text{ if } j\notin{S}_\lambda(\beta)\\ \widehat{\beta}_{j}^{\text{ols}} \text{ if } j\in{S}_\lambda(\beta),\end{array}\right.$$where $S_\lambda (\beta)$ refers to all non-zero coordinates when solving $(\ell0)$.

The problem $(\ell2)$ is strictly convex if $\ell$ is the quadratic norm, in other words, the Ridge estimator is always well defined, with in addition an explicit form for the estimator, $$\widehat{ {\beta}}_\lambda^{\text{ ridge}}=(\mathbf{X}^T\mathbf{X}+\lambda\mathbb{I})^{-1}\mathbf{X}^T\mathbf{y}=(\mathbf{X}^T\mathbf{X}+\lambda\mathbb{I})^{-1}(\mathbf{X}^T\mathbf{X})\widehat{ {\beta}}^{\text{ ols}}$$Therefore, it can be deduced that $$\text{bias}[\widehat{ {\beta}}_\lambda^{\text{ ridge}}]=-\lambda[\mathbf{X}^T\mathbf{X}+\lambda\mathbb{I}]^{-1}~\widehat{ {\beta}}^{\text{ ols}}$$and$$\text{Var}[\widehat{\beta}_\lambda^{\text{ ridge}}]=\sigma^2[\mathbf{X}^T\mathbf{X}+\lambda\mathbb{I}]^{-1}\mathbf{X}^T\mathbf{X}[\mathbf{X}^T\mathbf{X}+\lambda\mathbb{I}]^{-1}$$With a matrix of orthonormal explanatory variables (i.e. $\mathbf{X}^T \mathbf{X}=\mathbb{I}$), the expressions can be simplified $$\text{bias}[\widehat{ {\beta}}_\lambda^{\text{ ridge}}]=\frac{\lambda}{1+\lambda}~\widehat{ {\beta}}^{\text{ ols}}\text{ and }\text{Var}[\widehat{ {\beta}}_\lambda^{\text{ ridge}}]=\frac{\sigma^2}{(1+\lambda)^2}\mathbb{I}$$Observe that $\text{Var}[\widehat{ {\beta}}_\lambda^{\text{ ridge}}]<\text{Var}[\widehat{ {\beta}}^{\text{ ols}}]$. And because  $$\text{mse}[\widehat{ {\beta}}_\lambda^{\text{ ridge}}]=\frac{p\sigma^2}{(1+\lambda)^2}+\frac{\lambda^2}{(1+\lambda)^2}\beta^T\beta$$we obtain an optimal value for $\lambda$: $\lambda^\star=k\sigma^2/\beta^T\beta$

On the other hand, if $\ell$ is no longer the quadratic norm but the $\ell_1$ norm, the problem $(\ell1)$ is not always strictly convex, and in particular, the optimum is not always unique (for example if $\mathbf{X}^T \mathbf{X}$ is singular). But if it is strictly convex, then predictions $\mathbf{X}\beta$ will be unique. It should also be noted that two solutions are necessarily consistent in terms of sign of coefficients: it is not possible to have $\beta_j<0$ for one solution and $\beta_j>0$ for another. From a heuristic point of view, the program $(\ell1)$ is interesting because it allows to obtain in many cases a corner solution, which corresponds to a problem resolution of type $(\ell0)$ – as shown visually on Figure 2.

Figure 2 : Penalization based on norms $\ell_0$, $\ell_1$ and $\ell_2$ (from Hastie et al. (2016)).

Let us consider a very simple model: $y_i=x_i \beta+\varepsilon$, with a penalty $\ell_1$ and a loss function $\ell_2$. The problem $(\ell1)$ then becomes  $$\min\big\{\mathbf{y}^T\mathbf{y}-2\mathbf{y}^T\mathbf{x}\beta+\beta\mathbf{x}^T\mathbf{x}\beta+2\lambda|\beta|\big\}$$The first order condition is then $$-2\mathbf{y}^T\mathbf{x} + 2\mathbf{x}^T\mathbf{x}\widehat{\beta}\pm 2\lambda=0$$And the sign of the last term depends on the sign of $\beta$. Suppose that the least square estimator (obtained by setting $\lambda=0$) is (strictly) positive, i. e. $\mathbf{y}^T \mathbf{x}>0$. If $\lambda$ is not too big, we can imagine that $\beta$ is of the same sign as $\widehat{\beta}^{\text{mco}}$, and therefore the condition becomes $-2\mathbf{y}^T \mathbf{x}+2\mathbf{x}^T \mathbf{x}\beta+2\lambda=0$, and the solution is $$\widehat{\beta}_{\lambda}^{\text{ lasso}}=\frac{\mathbf{y}^T\mathbf{x}-\lambda}{\mathbf{x}^T\mathbf{x}}$$By increasing $\lambda$, there will be a time such that $\widehat{\beta}_λ=0$. If we increase $\lambda$ a bit little more, $\widehat{\beta}_λ$ does not become negative because in this case the last term of the first order condition changes, and in this case we try to solve $$-2\mathbf{y}^T\mathbf{x} + 2\mathbf{x}^T\mathbf{x}\widehat{\beta}- 2\lambda=0$$whose solution is then $$\widehat{\beta}_{\lambda}^{\text{ lasso}}=\frac{\mathbf{y}^T\mathbf{x}+\lambda}{\mathbf{x}^T\mathbf{x}}$$But this solution is positive (we assumed $\mathbf{y}^T \mathbf{x}>0$), and so it is possible to have $\widehat{\beta}_\lambda <0$at the same time. Also, after a while, $\widehat{\beta}_\lambda=0$, which is then a corner solution. Things are of course more complicated in larger dimensions (Tibshirani & Wasserman (2016) goes back at length on the geometry of the solutions) but as Candès & Plan (2009) notes, under minimal assumptions guaranteeing that the predictors are not strongly correlated, the Lasso obtains a quadratic error almost as good as if we had an oracle providing perfect information on the set of $\beta_j$‘s that are not zero. With some additional technical hypotheses, it can be shown that this estimator is “sparsistant” in the sense that the support of $\widehat{\beta}_\lambda^{\text{lasso}}$ is that of $\beta$, in other words Lasso has made it possible to select variables (more discussions on this point can be obtained in Hastie et al. (2016)).

More generally, it can be shown that $\widehat{\beta}_\lambda^{\text{lasso}}$ is a biased estimator, but may be of sufficiently low variance that the mean square error is lower than using least squares. To compare the three techniques, relative to the least square estimator (obtained when $\lambda=0$), if we assume that the explanatory variables are orthonormal, then $$\widehat{\beta}_{\lambda,j}^{\text{ subset}}=\widehat{\beta}_{j}^{\text{ ols}}\boldsymbol{1}_{|\widehat{\beta}_{\lambda,j}^{\text{ subset}}|>b}, ~~\widehat{\beta}_{\lambda,j}^{\text{ ridge}}=\frac{\widehat{\beta}_{j}^{\text{ ols}}}{1+\lambda}$$and$$\widehat{\beta}_{\lambda,j}^{\text{ lasso}}=\text{sign}[\widehat{\beta}_{j}^{\text{ ols}}]\cdot(|\widehat{\beta}_{j}^{\text{ ols}}|-\lambda)_+$$

Figure 3 : Penalization based on norms ,  and  (from Hastie et al. (2016)).

To be continued with probably a final post this week (references are online here)…

# Foundations of Machine Learning, part 3

This post is the seventh one of our series on the history and foundations of econometric and machine learning models. The first fours were on econometrics techniques. Part 6 is online here.

## Boosting and sequential learning

As we have seen before, modelling here is based on solving an optimization problem, and solving the problem described by equation $(6)$ is all the more complex because the functional space $\mathcal{M}$ is large. The idea of boosting, as introduced by Shapire & Freund (2012), is to learn, slowly, from the errors of the model, in an iterative way. In the first step, we estimate a model $m_1$ for $y$, from $\mathbf{X}$, which will give an error $\varepsilon_1$. In the second step, we estimate a model $m_2$ for $\varepsilon_1$, from $X$, which will give an error $\varepsilon_2$, etc. We will then retain as a model, after $k$ iterations $$m^{(k)}(\cdot)=\underbrace{m_1(\cdot)}_{\sim y}+\underbrace{m_2(\cdot)}_{\sim \epsilon_1}+\underbrace{m_3(\cdot)}_{\sim \epsilon_2}+\cdots+\underbrace{m_k(\cdot)}_{\sim \epsilon_{k-1}}=m^{(k-1)}(\cdot)+m_k(\cdot)~~~(7)$$Here, the error $\varepsilon$ is seen as the difference between $y$ and the model $m(\mathbf{x})$, but it can also be seen as the gradient associated with the quadratic loss function. Formally, $\varepsilon$ can be seen as $\nabla\ell$ in a more general context (here we find an interpretation that reminds us of residuals in generalized linear models).

Equation $(7)$ can be seen as a descent of the gradient, but written in a dual way. The problem will then be rewritten as an optimization problem: $$m^{(k)}=m^{(k-1)}+\underset{h\in\mathcal{H}}{\text{argmin}}\left\lbrace \sum_{i=1}^n \ell(\underbrace{y_i-m^{(k-1)}(\boldsymbol{x}_i)}_{\varepsilon_{k,i}},h(\boldsymbol{x}_i))\right\rbrace~~~(8)$$where the trick is to consider a relatively simple space $\mathcal{H}$ (we will speak of “weak learner”). Classically, $\mathcal{H}$ functions are step-functions (which will be found in classification and regression trees) called “stumps”. To ensure that learning is indeed slow, it is not uncommon to use a shrinkage parameter, and instead of setting, for example, $\varepsilon_1=y-m_1 (\mathbf{x})$, we will set $\varepsilon_1=y-\alpha\cdot m_1 (\mathbf{x})$ with $\alpha\in[0.1]$. It should be noted that it is because a non-linear space is used for $\mathcal{H}$, and learning is slow, that this algorithm works well. In the case of the Gaussian linear model, remember that the residuals $\varepsilon=y-\mathbf{x}^T\beta$ are orthogonal to the explanatory variables, $\mathbf{X}$, and it is then impossible to learn from our errors. The main difficulty is to stop in time, because after too many iterations, it is no longer the m function that is approximated, but the noise. This problem is called overlearning.

This presentation has the advantage of having a heuristic reminiscent of an econometric model, by iteratively modelling the residuals by a (very) simple model. But this is often not the presentation used in the learning literature, which places more emphasis on an optimization algorithm heuristic (and gradient approximation). The function is learned iteratively, starting from a constant value, $$m^{(0)}=\underset{m\in\mathbb{R}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \ell(y_i,m)\right\rbrace$$then we consider the following learning procedure$${\displaystyle m^{(k)}=m^{(k-1)}+{\underset{h\in {\mathcal {H}}}{\text{argmin}}}\sum _{i=1}^{n}\ell(y_{i},m^{(k-1)}(\mathbf{x}_{i})+h(\mathbf{x}_{i}))}~~~(9)$$which can be written, if $\mathcal{H}$ is a set of differentiable functions,$${\displaystyle m^{(k)}=m^{(k-1)}-\gamma_{k}\sum _{i=1}^{n}\nabla _{m^{(k-1)}}\ell(y_{i},m^{(k-1)}(\mathbf{x}_{i})),}$$where $${\displaystyle \gamma _{k}=\underset{\gamma }{\text{argmin }}\sum _{i=1}^{n}\ell\left(y_{i},m^{(k-1)}( \mathbf{x}_{i})-\gamma \nabla _{m^{(k-1)}}\ell(y_{i},m^{(k-1)}( \mathbf{x}_{i}))\right).}$$To better understand the relationship with the approach described above, at step $k$, pseudo-residuals are defined by setting $$r_{i,k}=-\left.\frac{\partial \ell(y_i,m(\mathbf{x}_i))}{\partial m(\mathbf{x}_i)}\right\vert_{m(\mathbf{x})=m^{(k-1)}( \mathbf{x})}\text{ where }i=1,\cdots,n$$A simple model is then sought to explain these pseudo-residuals according to the explanatory variables $\mathbf{x}_i$, i.e. $r_{i,k}=h^\star(\mathbf{x}_i)$, where $h^\star\in\mathcal{H}$. In a second step, we look for an optimal multiplier by solving$$\gamma_k = \underset{\gamma\in\mathbb{R}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \ell(y_i,m^{(k-1)}( \mathbf{x}_i)+\gamma h^\star(\mathbf{x}_i))\right\rbrace$$then update the model by setting $$m_k (\cdot)=m_(k-1) (\cdot)+\gamma_k h^\star (\cdot)$$. More formally, we move from equation $(8)$ – which clearly shows that we are building a model on residuals – to equation $(9)$ – which will then be translated as a gradient calculation problem – noting that $\ell(y,m+h)=\ell(y-m,h)$. Classically, class $\mathcal{H}$ of functions consists in regression trees. It is also possible to use a form of penalty by setting $m_k (\cdot)=m_(k-1) (\cdot)+\nu\gamma_k h^\star (\cdot)$, with $\nu\in(0,1)$. But let’s go back  a little further – in our next post – on the importance of penalization before discussing the numerical aspects of optimization.

To be continued (keep in mind that references are online here)…

# Foundations of Machine Learning, part 2

This post is the sixth one of our series on the history and foundations of econometric and machine learning models. The first fours were on econometrics techniques. Part 5 is online here.

## The probabilistic formalism in the 80’s

We have a training sample, with observations $(\mathbf{x}_i,y_i)$ where the variables $y$ are in a set $\mathcal{Y}$. In the case of classification, $\mathcal{Y}=\{-1,+1\}$, but a relatively general set can be considered (note that if econometricians prefer $\mathcal{Y}=\{0,1\}$ – because of the Bernoulli distribution and because $0$ and $1$ are lower and upper bounds of probabilities, people in the “machine learning” community prefer $\mathcal{Y}=\{-1,+1\}$). A predictor $m$ is an function taking values in $\mathcal{Y}$, used to label (or classify) future new observations, using some features that lie in a set $\mathcal{X}$. It is assumed that the labels are produced by an (unknown) classifier $f$ called target. For a statistician, this function would be the real model. Naturally, we want to build $m$ as close as possible to $f$. Let $\mathbb{P}$ be a (unknown) distribution on $\mathcal{X}$. The error of $m$ with respect to target $f$ is defined by $$\mathcal{R}_{\mathbb{P},f}(m)=\mathbb{P}[m(\boldsymbol{X})\neq f(\boldsymbol{X})]\text{ where }\boldsymbol{X}\sim\mathbb{P}$$or equivalently,$$\mathcal{R}_{\mathbb{P},f}(m)=\mathbb{P}\big[\{\boldsymbol{x}\in\mathcal{X}:m(\boldsymbol{x})\neq f(\boldsymbol{x})\}\big]$$To obtain our “optimal” classifier, it becomes necessary to assume that there is a link between the data in our sample and the pair $(\mathbb{P},f)$, i.e. a data generation model. We will then assume that the $\mathbf{x}_i$ are obtained by independent draws according to $\mathbb{P}$, and that then $y_i=f(\mathbf{x}_i)$ . We can define the empirical risk of a classifier $m$, as $$\widehat{{R}}(m)=\frac{1}{n}\sum_{i=1}^n \boldsymbol{1}(m(\boldsymbol{x}_i)\neq y_i)$$

It is important to recognize that a perfect model cannot be found, in the sense that $R_{\mathbb{P},f} (m)=0$. Indeed, if we consider the simplest case, with $\mathcal{X}=\{x_1,x_2\}$ and $\mathbb{P}$ is such that $\mathbb{P}(\{x_1\})=p$ and $\mathbb{P}(\{x_2\})=1-p$. The probability of never observing $\{x_2\}$ among the $n$ observations is $(1-p)^n$, and if $p<1/n$, it is quite likely never to observe $\{x_2\}$ so it can never be predicted. We cannot therefore hope to have a zero risk whatever $\mathbb{P}$. And more generally, it is also possible to observe $\{x_1\}$ and $\{x_2\}$, and despite everything, to make mistakes on the labels. Also, instead of looking for a perfect model, we can try to have an “approximately correct” model. We will then try to find $m$ such that $R_{\mathbb{P},f} (m)\leq\varepsilon$, where $\varepsilon$ is an a priori specified threshold. But even this condition is too strong, and cannot be fulfilled. Thus, we will usually as to have $R_{\mathbb{P},f} (m)\leq\varepsilon$ with some probability $1-\delta$. Hence, we will try to be “probably approximately correct” (PAC), allowing to make a mistake with a probability $\delta$, again fixed a priori.

Also, when we build a classifier, we do not know either $\mathbb{P}$ or $f$, but we give ourselves a precision criterion $\varepsilon$, and a confidence parameter $\delta$, and we have $n$ observations. Note that $n$, $\varepsilon$ and $\delta$ can be linked. We then look for a model $m$ such that $R_{\mathbb{P},f} (m)\leq\varepsilon$ with probability (at least) $1-\delta$, so that we are probably approximately correct. Wolpert (1996) has shown (see details in Wolpert & Macready (1997)) that there is no universal learning algorithm. In particular, it can be shown that there is $\mathbb{P}$ such that $R_{\mathbb{P},f} (m)$ is relatively high, with a relatively high probability (also).

The interpretation is that since we cannot learn (in the PAC sense) about all the functions $m$, we will then force $m$ to belong to a particular class, noted $\mathcal{M}$. Let us suppose, to start with, that $\mathcal{M}$ contains a finite number of possible models. We can then show that for all $\varepsilon$ and $\delta$, that for all $\mathbb{P}$ and $f$, if we have enough observations (more precisely $n\geq \varepsilon^{-1} \log[\delta^{-1} |\mathcal{M}|]$, then with a greater probability than $1-\delta$, $R_{\mathbb{P},f} (m^\star)\leq\varepsilon$ where$$m^\star \in \underset{m\in\mathcal{M}}{\text{argmin}}\Big\lbrace\frac{1}{n}\sum_{i=1}^n \boldsymbol{1}(m(\boldsymbol{x}_i)\neq y_i)\Big\rbrace$$in other words $m^\star$ is a model in $\mathcal{M}$ that minimizes empirical risk.

We can go a little further, staying in the case where $\mathcal{Y}=\{-1,+1\}$. An $\mathcal{M}$ class of classifiers will be called PAC-learnable if there is $n_M:[0,1]^2\rightarrow \mathbb{N}$ such that, for all $\varepsilon$, $\delta$, $\mathbb{P}$ and if it is assumed that the target $f$ belongs to $\mathcal{M}$, then using $n>n_M (\varepsilon,\delta)$ observations $\mathbf{x}_i$ drawn from $\mathbb{P}$, labelled $y_i$ by $f$, then there is $m\in\mathcal{M}$ such that, with probability $1-\delta$, $R_{\mathbb{P},f} (m)\leq\varepsilon$. The $n_M$ function is then called “sample complexity to learn”. In particular, we have seen that if $M$ contains a finite number of classifiers, then $\mathcal{M}$ is PAC-learnable with complexity $n_M (\varepsilon,\delta)=\varepsilon^{-1} \log[\delta^{-1} |M|]$.

Naturally, we would like to have a more general result, especially if $\mathcal{M}$ is not finite. To do this, the $VC$ dimension of Vapnik-Chervonenkis must be used, which is based on the idea of shattering points (for a binary classification). Consider $k$ points $\{x_1,\cdot,x_k\}$, and consider the set $${E}_k=\big\lbrace(m(\boldsymbol{x}_1),\cdots,m(\boldsymbol{x}_k))\text{ for }m\in\mathcal{M})\big\rbrace$$ Note that the elements of $E_k$ belong to $\{-1,+1\}^k$. In other words, $|E_k |\leq 2^k$. We will say that M shatter all the points if all the combinations are possible, i. e. $|E_k |=2^k$. Intuitively, the labels of the set of points do not provide enough information on target $f$, because anything is possible. The $VC$ dimension of $\mathcal{M}$ is then$$VC(\mathcal{M})=\sup\big\lbrace k\text{ such that }\mathcal{M}\text{ shatters }\{\boldsymbol{x}_1,\cdots\boldsymbol{x}_k\}\big\rbrace$$

For example, if $\mathcal{X}=\mathbb{R}$ and all (simple) models of the form [1] $m_{a,b}=\mathbf{1}_{\pm}(x\in[a,b])$ are considered. No set of $\{x_1,x_2,x_2,x_3\}$ ordered points can be shattered because it is sufficient to assign respectively +1, -1 and +1 to $x_1$, $x_2$ and $x_3$ respectively, therefore $VC<3$. On the other hand $\{0,1\}$ is shattered, so $VC\geq 2$. The dimension of this predictor set is $2$: If we increase by one dimension, $\mathcal{X}=\mathbb{R}^2$ and consider all (simple) models of the form $m_{a,b}=\mathbf{1}_{\pm} (x\in[a,b])$ (where $[a,b]$ refers to the rectangle), then the dimension of $\mathcal{M}$ is here $4$.

To introduce SVMs, let’s place ourselves in the case where $\mathcal{X}=\mathbb{R}^k$, and consider separations by hyperplanes passing through the origin (we will say homogeneous), in the sense that $m_{\mathbf{w}} (\mathbf{x})=\mathbf{1}_{\pm}(\mathbf{w}^T \mathbf{x}\geq 0)$. It can be shown that no set of $k+1$ points can be shattered by these two homogeneous spaces in $\mathbb{R}^k$, and therefore $VC(M)=k$. If we add a constant, in the sense that $m_{\mathbf{w},b} (\mathbf{x})=\mathbf{1}_{\pm}(\mathbf{w}^T \mathbf{x}+b\geq 0)$, we can show that no set of $k+2$ points can be sprayed by these two (non-homogeneous) spaces in $\mathbb{R}^k$, and therefore $VC(M)=k+1$. This dimension reminds us of the dimension of the model we’ve seen in the econometric context.

From this dimension $VC$, we deduce the so-called fundamental theorem of learning: if $\mathcal{M}$ is a class of dimension $d=VC(M)$, then there are positive constants $\underline{C}$ and $\overline{C}$ such as the sample complexity for M to be PAC-learnable satisfies$$\underline{C}\epsilon^{-1}\big(d+\log[\delta^{-1}]\big)\leq n_{\mathcal{M}}(\epsilon,\delta) \leq \overline{C}\epsilon^{-1}\big(d\log[\epsilon^{-1}]+\log[\delta^{-1}]\big)$$The link between the notion of learning (as defined in Vailiant (1984)) and the $VC$ dimension was clearly established in Blumer et al (1989).

Nevertheless, while the work of Vapnik and Chervonenkis is considered to be the foundation of statistical learning, Thomas Cover’s work in the 1960s and 1970s should also be mentioned, in particular Cover (1965) on the capacities of linear models, and Cover & Hart (1967) on learning in the context of the algorithm of the $k$-nearest neighbors. These studies have linked learning, information theory (with the textbook Cover & Thomas (1991)), complexity and statistics. Other authors have subsequently brought the two communities closer together, in terms of learning and statistics. For example, Halbert White proposed to see neural networks in a statistical context in White (1989), going so far as to state that « learning procedures used to train artificial neural networks are inherently statistical techniques. It follows that statistical theory can provide considerable insight into the properties, advantages, and disadvantages of different network learning methods ». This turning point in the late 1980s will anchor learning theory in a probabilistic context.

## Objective and loss function

These choices (of objective and loss function) are essential, and very dependent on the problem under consideration. Let us begin by describing a historically important model, Rosenblatt’s (1958) “perceptron”, introduced into classification problems, where $y\in\{-1,+1\}$, inspired by McCulloch & Pitts (1943). We have data $\{(y_i,\mathbf{x}_i)\}$, and we will iteratively build a set of $m_k[\mathbf{x}$ models, where at each step, we will learn from the errors of the previous model. In the perceptron, a linear model is considered so that :$$m(\mathbf{x})=\boldsymbol{1}_{\pm}(\beta_0+\mathbf{x}^T \boldsymbol{\beta}\geq 0)=\left\lbrace\begin{array}{l}+1\text{ si }\beta_0+\mathbf{x}^T \boldsymbol{\beta}\geq 0\\-1\text{ si }\beta_0+\mathbf{x}^T \boldsymbol{\beta}< 0\end{array}\right.$$where $\beta$ coefficients are often interpreted as “weights” assigned to each of the explanatory variables. We give ourselves initial weights $(\beta_0^{(0)},\beta^{(0)}$, which we will update taking into account the prediction error made, between $y_i$ and the prediction $\widehat{y}_i^{(k)}$ :$$\widehat{y}_i^{(k)}=m^{(k)}(\mathbf{x}_i)=\boldsymbol{1}_{\pm}(\beta_0^{(k)}+\mathbf{x}^T \boldsymbol{\beta}^{(k)}\geq 0),$$with, in the case of the perceptron:$$\beta_j^{(k+1)}={\beta}_j^{(k)}+\eta\underbrace{(\mathbf{y}-\widehat{\mathbf{y}}^{(k)})^T}_{=\ell({\mathbf{y}},\widehat{\mathbf{y}}^{(k)})}\mathbf{x}_j$$Here $\ell(y,y')=\mathbf{1}(y\neq y')$ is a loss function, which will allow to give a price to an error made, by predicting $\widehat{y}=m(\mathbf{x})$ and observing $y$. For a regression problem, we can consider a quadratic error $\ell_2$, such that $\ell(y,m(\mathbf{x}))=(y-m(\mathbf{x}))^2$ or in absolute value $\ell_1$, with $\ell(y,m(\mathbf{x}))=|y-m(\mathbf{x})|$. Here, for our classification problem, we used a mis-qualification indicator (we could discuss the symmetry of this loss function, suggesting that a false positive costs as much as a false negative). Once this loss function has been specified, we recognize in the problem previously described a gradient descent, and we see that we are trying to solve:$$m^\star(\mathbf{x})=\underset{m\in\mathcal{M}}{\text{argmin}}\left\lbrace\sum_{i=1}^n \ell(y_i,m(\mathbf{x}_i))\right\rbrace~~~(6)$$for a predefined set of predictors $\mathcal{M}$. Any machine learning problem is mathematically formulated as an optimization problem, whose solution determines a set of model parameters (if the $\mathcal{M}$ family is described by a set of parameters – which can be coordinates in a functional database). We can note $\mathcal{M}_0$ the space of the hyperplanes of $\mathbb{R}^p$ in the sense that$$m\in\mathcal{M}_0 \text{\quad means \quad}m(\mathbf{x})=\beta_0+\beta^T\mathbf{x}\text{ where }\beta\in\mathbb{R}^p$$ generating the class of linear predictors. We will then have the estimator that minimizes the empirical risk. Some of the recent work in statistical learning aims to study the properties of the estimator $\widehat{m}^\star$, known as “oracle”, in a family of $\mathcal{M}$ estimators, $$\widehat{m}^{\star} =\underset{\widehat{m}\in\mathcal{M}}{\text{argmin}}\big\lbrace\mathcal{R}(\widehat{m},m)\big\rbrace$$This estimator is, of course, impossible to define because it depends on $m$, the real model, unknown.

But let’s come back a little more to these loss functions. A loss function $\ell$ is a function $\mathbb{R}^d\times\mathbb{R}^d\rightarrow\mathbb{R}_+$, symmetric, which checks the triangular inequality, and such that $\ell(x,y)=0$ if and only if $x=y$. The associated norm is $\|\cdot\|$, such that $\ell(x,y)=\|x-y\|=\ell(x-y,0)$ (using the fact that $\ell(x,y+z)=\ell(x-y,z)$ – we will review this fundamental property later).

For a quadratic loss function, it should be noted that we can have a particular interpretation of this problem, since:$$\overline{y}=\underset{m\in\mathbb{R}}{\text{argmin}} \left\lbrace\sum_{i=1}^n\frac{1}{n} [y_i-m]^2\right\rbrace=\underset{m\in\mathbb{R}}{\text{argmin}} \left\lbrace \sum_{i=1}^n \ell_2(y_i,m)\right\rbrace$$ where $\ell_2$ is the usual quadratic distance If we assume – as we did in econometrics – that there is an underlying probabilistic model, and observe that : $$\displaystyle{\mathbb{E}(Y)=\underset{m\in\mathbb{R}}{\text{argmin}}\left\lbrace\mathbb{E}\left([Y-m]^2\right)\right\rbrace=\underset{m\in\mathbb{R}}{\text{argmin}}\left\lbrace\mathbb{E}\big[\ell_2(Y,m)\big]\right\rbrace}$$it should be noted that what we are trying to obtain here, by solving the problem $(6)$ by taking the norm $\ell_2$, is an approximation (in a given functional space, $\mathcal{M}$) of the conditional expectation $x\mapsto\mathbb{E}[Y|\mathbf{X}=\mathbf{x}]$. Another particularly interesting loss function is the loss $\ell_1$,$\ell_1 (y,m)=|y-m|[\latex]. It should be recalled that [latex display="true"]\displaystyle{\text{median}(\boldsymbol{y})=\underset{m\in\mathbb{R}}{\text{argmin}}\left\lbrace\sum_{i=1}^n\ell_1(y_i,m)\right\rbrace}$The optimization problem :$$\widehat{m}^{\star}=\underset{m\in\mathcal{M}_0}{\text{argmin}}\left\lbrace\sum_{i=1}^n\vert y_i-m(\mathbf{x}_i)\vert\right\rbrace$$ is obtained in econometrics by assuming that the conditional law of $Y$ follows a Laplace law centered on $m(\mathbf{x})$, and by maximizing the likelihood (log) (the sum of the absolute values of the errors corresponds to the log-reasonableness of a Laplace law). It should also be noted that if the conditional law of $Y$ is symmetrical with respect to $0$, the median and the mean coincide If this loss function is rewritten  $$\ell_1(y,m)=\vert (y-m)(1/2-\boldsymbol{1}_{y\leq m})\vert$$ a generalization can be obtained for $\tau\in[0.1]$:$$\widehat{m}^\star_\tau=\underset{m\in\mathcal{M}_0}{\text{argmin}}\left\lbrace\sum_{i=1}^n \ell_\tau^{ q} (y_i,m(\mathbf{x}_i)) \right\rbrace$$where$$\ell_{\tau}^{q}(x,y)= (x-y)(\tau-\boldsymbol{1}_{x\leq y})$$  is then the quantile regression of level $\tau$ (Koenker, 2003; d'Haultefœuille & Givord, 2014). Another loss function, introduced by Aigner et al (1977) and analysed in Waltrup et al (2014), is the function associated with the notion of expectations: $$\displaystyle{\ell}^{\text{ e}}_{\tau}(x,y)= (x-y)^2\cdot\big\vert\tau-\boldsymbol{1}_{x\leq y}\big\vert$$with $\tau\in[0.1]$. We see the parallel with the quantile function: $$\displaystyle{\ell}^{\text{ q}}_{\tau}(x,y)= \vert x-y\vert \cdot\big\vert\tau-\boldsymbol{1}_{x\leq y}\big\vert$$Koenker & Machado (1999) and Yu & Moyeed (2001) also noted a link between this condition and the search for maximum likelihood when $Y$'s conditional law follows an asymmetric Laplace law.

In connection with this approach, Gneiting (2011) introduced the notion of "ellicable statistics" - or "ellicable measurement" in its probabilistic (or distributional) version: a statistic $T$ will be said to be "ellicitable" if there is a loss function $\ell:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}_+$ such that:$$T(Y)=\underset{x\in\mathbb{R}}{\text{argmin}}\left\lbrace\int_{\mathbb{R}} \ell(x,y)dF(y)\right\rbrace=\underset{x\in\mathbb{R}}{\text{argmin}}\left\lbrace\mathbb{E}\big[ \ell(x,Y)\big]\text{ where }Y\overset{\mathcal{L}}{\sim} F\right\rbrace$$ The mean (mathematical expectation) is thus ellicable by the quadratic distance, $\ell_2$, while the median is ellicable by the distance $\ell_1$. According to Gneiting (2011), this property is essential for obtain predictions and forecasts. There may then be a strong link between measures associated with probabilistic models and loss functions. Finally, Bayesian statistics provide a direct link between the form of the a priori law and the loss function, as studied by Berger (1985) and Bernardo & Smith (2000). We will come back to the use of these different norms in the section on penalization.

To be continued (keep in mind that references are online here)…

[1] Where the indicator $\mathbf{1}_{\pm}$ does not take values 0 or 1 (like the classical $\mathbf{1}$ function), but -1 and +1.

# References on Econometrics and Machine Learning

In our series of posts on the history and foundations of econometric and machine learning models, a lot of references where given. Here they are.

# Foundations of Machine Learning, part 1

This post is the fifth one of our series on the history and foundations of econometric and machine learning models. The first fours were on econometrics techniques. Part 4 is online here.

In parallel with these tools developed by, and for economists, a whole literature has been developed on similar issues, centered on the problems of prediction and forecasting. For Breiman (2001a), a first difference comes from the fact that the statistic has developed around the principle of inference (or to explain the relationship linking $y$ to variables $\mathbf{x}$) while another culture is primarily interested in prediction. In a discussion that follows the article, David Cox states very clearly that in statistic (and econometrics) “predictive success (…) is not the primary basis for model choice“. We will get back here on the roots of automatic learning techniques. The important point, as we will see, is that the main concern of machine learning is related to the generalization properties of a model, i.e. its performance – according to a criterion chosen a priori – on new data, and therefore on non-sample tests.

## A learning machine

Today, we speak of “machine learning” to describe a whole set of techniques, often computational, as alternatives to the classical econometric approach. Before characterizing them as much as possible, it should be noted that historically other names have been given. For example, Friedman (1997) proposes to make the link between statistics (which closely resemble econometric techniques – hypothesis testing, ANOVA, linear regression, logistics, GLM, etc.) and what was then called “data mining” (which then included decision trees, methods from the closest neighbours, neural networks, etc.). The bridge between those two cultures corresponds to “statistical learning” techniques described in Hastie et al (2009). But one should keep in mind that machine learning is a very large field of research.

The so-called “natural” learning (as opposed to machine learning) is that of children, who learn to speak, read and play. Learning to speak means segmenting and categorizing sounds, and associating them with meanings. A child also learns simultaneously the structure of his or her mother tongue and acquires a set of words describing the world around him or her. Several techniques are possible, ranging from rote learning, generalization, discovery, more or less supervised or autonomous learning, etc. The idea in artificial intelligence is to take inspiration from the functioning of the brain to learn, to allow “artificial” or “automatic” learning, by a machine. A first application was to teach a machine to play a game (tic-tac-toe, chess, go, etc.). An essential step is to explain the objective it must achieve to win. One historical approach has been to teach the machine the rules of the game. If it allows you to play, it will not help the machine to play well. Assuming that the machine knows the rules of the game, and that it has a choice between several dozen possible moves, which one should it choose? The classical approach in artificial intelligence uses the so-called min-max algorithm using an evaluation function: in this algorithm, the machine searches forward in the possible moves tree, as far as the calculation resources allow (about ten moves in chess, for example). Then, it calculates different criteria (which have been previously indicated to her) for all positions (number of pieces taken, or lost, occupancy of the center, etc. in our example of the chess game), and finally, the machine plays the move that allows it to maximize its gain. Another example may be the classification and recognition of images or shapes. For example, the machine must identify a number in a handwritten handwriting (checks, ZIP code on envelopes, etc). It is a question of predicting the value of a variable $y$, knowing that a priori $y\in\{0,1,2,\cdots,8,9\}$. A classical strategy is to provide the machine with learning bases, in other words here millions of labelled (identified) images of handwritten numbers. A simple (and natural) strategy is to use a decision criterion based on the closest neighbors whose labels are known (using a predefined metric).

The method of the closest neighbors (“$k$-nearest neighbors”) can be described as follows: we consider (as in the previous part) a set of n observations, i. e. pairs $(y_i,\mathbf{x}_i)$ with $\mathbf{x}_i\in\mathbb{R}^p$. Let us consider a distance $\Delta$ on $\mathbb{R}^p$ (the Euclidean distance or the Mahalanobis distance, for example). Given a new observation $\mathbf{x}\in\mathbb{R}^p$, let us assume the ordered observations as a function of the distance between the $\mathbf{x}_i$ and $\mathbf{x}$, in the sense that $$\Delta(\mathbf{x}_1, \mathbf{x})\leq\Delta(\mathbf{x}_2, \mathbf{x})\leq\cdots\leq\Delta(\mathbf{x}_n, \mathbf{x})$$ then we can consider as prediction for y the average of the nearest $k$ neighbours,$$\widehat{m}_k(\mathbf{x})=\frac{1}{k}\sum_{i=1}^k y_i$$Learning here works by induction, based on a sample (called the learning – or training – sample).

Automatic learning includes those algorithms that give computers the ability to learn without being explicitly programmed (as Arthur Samuel defined it in 1959). The machine will then explore the data with a specific objective (such as searching for the nearest neighbours in the example just described). Tom Mitchell proposed a more precise definition in 1998: a computer program is said to learn from experience $E$ in relation to a task $T$ and a performance measure $P$, if its performance on $T$, measured by $P$, improves with experience $E$. Task $T$ can be a defect score for example, and performance $P$ can be the percentage of errors made. The system learns if the percentage of predicted defects increases with experience.

As we can see, machine learning is basically a problem of optimizing a criterion based on data (from now on called learning). Many textbooks on machine learning techniques propose algorithms, without ever mentioning any probabilistic model. In Watt et al (2016) for example, the word “probability” is mentioned only once, with this footnote that will surprise and make smile any econometricians, “the logistic regression can also be interpreted from a probabilistic perspective” (page 86). But many recent books offer a review of machine learning approaches using probabilistic theories, following the work of Vaillant and Vapnik. By proposing the paradigm of “probably almost correct” learning (PAC), a probabilistic flavor has been added to the previously very computational approach, by quantifying the error of the learning algorithm (usually in a classification problem).

To be continued (references are online here)…

# Probabilistic Foundations of Econometrics, part 4

This post is the fourth one of our series on the history and foundations of econometric and machine learning models. Part 3 is online here.

## Goodness of Fit, and Model

In the Gaussian linear model, the determination coefficient – noted $R^2$ – is often used as a measure of fit quality. It is based on the variance decomposition formula $$\underbrace{\frac{1}{n}\sum_{i=1}^n (y_i-\bar{y})^2}_{\text{total variance}}=\underbrace{\frac{1}{n}\sum_{i=1}^n (y_i-\widehat{y}_i)^2}_{\text{residual variance}}+\underbrace{\frac{1}{n}\sum_{i=1}^n (\widehat{y}_i-\bar{y})^2}_{\text{explained variance}}$$ The $R^2$ is defined as the ratio of explained variance and total variance, another interpretation of the coefficient that we had introduced from the geometry of the least squares $$R^2= \frac{\sum_{i=1}^n (y_i-\bar{y})^2-\sum_{i=1}^n (y_i-\widehat{y}_i)^2}{\sum_{i=1}^n (y_i-\bar{y})^2}$$The sums of the error squares in this writing can be rewritten as a log-likelihood. However, it should be remembered that, up to one additive constant (obtained with a saturated model) in generalized linear models, deviance is defined by ${Deviance}(\widehat{\beta}) = -2\log[\mathcal{L}]$ which can also be noted $Deviance(\widehat{\mathbf{y}})$. A null deviance can be defined as the one obtained without using the explanatory variables $\mathbf{x}$, so that $\widehat{y}_i=\overline{y}$. It is then possible to define, in a more general context (with a non-Gaussian distribution for $y$)$$R^2=\frac{{Deviance}(\overline{y})-{Deviance}(\widehat{\mathbf{y}})}{{Deviance}(\overline{y})}=1-\frac{{Deviance}(\widehat{\mathbf{y}})}{{Deviance}(\overline{y})}$$However, this measure cannot be used to choose a model, if one wishes to have a relatively simple model in the end, because it increases artificially with the addition of explanatory variables without significant effect. We will then tend to prefer the adjusted $R^2$,$$\bar R^2 = {1-(1-R^{2})\cdot{n-1 \over n-p}} = R^{2}-\underbrace{(1-R^{2})\cdot{p-1 \over n-p}}_{\text{penalty}}$$where $p$ is the number of parameters of the model. Measuring the quality of fit will penalize overly complex models.

This idea will be found in the Akaike criterion, where $AIC=Deviance+2\cdot p$ or in the Schwarz criterion, $BIC=Deviance+log(n)\cdot p$. In large dimensions (typically $p>\sqrt{n}$), we will tend to use a corrected $AIC$, defined by $AIC_c=Deviance+2⋅p⋅n/(n-p-1)$.

These criterias are used in so-called “stepwise” methods, introducing the set methods. In the “forward” method, we start by regressing to the constant, then we add one variable at a time, retaining the one that lowers the $AIC$ criterion the most, until adding a variable increases the AIC criterion of the model. In the “backward” method, we start by regressing on all variables, then we remove one variable at a time, removing the one that lowers the $AIC$ criterion the most, until removing a variable increases the $AIC$ criterion from the model.

Another justification for this notion of penalty (we will come back to this idea in machine learning) can be the following. Let us consider an estimator in the class of linear predictors,$$\mathcal{M}=\big\lbrace m:~m(\mathbf{x})=s_h(\mathbf{x})^T\mathbf{y} \text{ where }S=(s(\mathbf{x}_1),\cdots,s(\mathbf{x}_n))^T\text{ is some smoothing matrix}\big\rbrace$$ and assume that $y=m_0 (x)+\varepsilon$, with $\mathbb{E}[\varepsilon]=0$ and $Var[\varepsilon]=\sigma^2\mathbb{I}$, so that $m_0 (x)=\mathbb{E}[Y|X=x]$. From a theoretical point of view, the quadratic risk, associated with an estimated model $\widehat{m}$, $\mathbb{E}\big[(Y-\widehat{m}(\mathbf{X}))^2\big]$, is written$$\mathcal{R}(\widehat{m})=\underbrace{\mathbb{E}\big[(Y-m_0(\mathbf{X}))^2\big]}_{\text{error}}+\underbrace{\mathbb{E}\big[(m_0(\mathbf {X})-\mathbb{E}[\widehat{m}(\mathbf{X})])^2\big]}_{\text{bias}^2}+\underbrace{\mathbb{E}\big[(\mathbb{E}[\widehat{m}(\mathbf{X})]-\widehat{m}(\mathbf{X}))^2\big]}_{\text{variance}}$$if $m_0$ is the true model. The first term is sometimes called “Bayes error”, and does not depend on the estimator selected, $\widehat{m}$.

The empirical quadratic risk, associated with a model $m$, is here:$$\widehat{\mathcal{R}}_n(m)=\frac{1}{n}\sum_{i=1}^n (y_i-m(\mathbf{x}_i))^2$$ (by convention). We recognize here the mean square error, “mse”, which will more generally give the “risk” of the model $m$ when using another loss function (as we will discuss later on). It should be noted that:$$\displaystyle{\mathbb{E}[\widehat{\mathcal{R}}_n(m)]=\frac{1}{n}\|m_0(\mathbf{x})-m(\mathbf{x})\|^2+\frac{1}{n}\mathbb{E}\big(\|{Y}-m_0(\mathbf{X})\|^2\big)}$$We can show that:$$n\mathbb{E}\big[\widehat{\mathcal{R}}_n(\widehat{m})\big]=\mathbb{E}\big(\|Y-\widehat{m}(\mathbf{x})\|^2\big)=\|(\mathbb{I}-\mathbf{S})m_0\|^2+\sigma^2\|\mathbb{I}-\mathbf{S}\|^2$$so that the (real) risk of $\widehat{m}$ is:$${\mathcal{R}}_n(\widehat{m})=\mathbb{E}\big[\widehat{\mathcal{R}}_n(\widehat{m})\big]+2\frac{\sigma^2}{n}\text{trace}(\boldsymbol{S})$$So, if $\text{trace}(\boldsymbol{S})\geq0$ (which is not a too strong assumption), the empirical risk underestimates the true risk of the estimator. Actually, we recognize here the number of degrees of freedom of the model, the right-hand term corresponding to Mallow’s $C_p$, introduced in Mallows (1973) using not deviance but $R^2$.

## Statistical Tests

The most traditional test in econometrics is probably the significance test, corresponding to the nullity of a coefficient in a linear regression model. Formally, it is the test of $H_0:\beta_k=0$ against $H_1:\beta_k\neq 0$. The so-called Student test, based on the statistics $t_k=\widehat{\beta}_k/se_{\widehat{β}_k}$, allows to decide between the two alternatives, using the test $p$-value, defined by $\mathbb{P}[|T|>|t_k|]$ avec $T\overset{\mathcal{L}}{\sim} Std_\nu$, where $\nu$ is the number of degrees of freedom of the model ($\nu=p+1$ for the standard linear model). In large dimension, however, this statistic is of very limited interest, given a significant FDR (“False Discovery Ratio”). Classically, with a level of significance $\alpha=0.05$, 5% of the variables are falsely significant. Suppose that we have $p$=100 explanatory variables, but that 5 (only) are really significant. We can hope that these 5 variables will pass the Student test, but we can also expect that 5 additional variables (false positive test) will emerge. We will then have 10 variables perceived as significant, while only half are significant, i.e. an FDR ratio of 50%. In order to avoid this recurrent pitfall in multiple tests, it is natural to use the procedure of Benjamini & Hochberg (1995).

## From a correlation to some causal effect

Econometric models are used to implement public policy evaluations. It is therefore essential to fully understand the underlying mechanisms in order to know which variables actually make it possible to act on a variable of interest. But then we move on to another important dimension of econometrics. Jerry Neyman was responsible for the first work on the identification of causal mechanisms, and then Rubin (1974) formalized the test, called the “Rubin causal model” in Holland (1986). The first approaches to the notion of causality in econometrics were based on the use of instrumental variables, models with discontinuity of regression, analysis of differences in differences, and natural or unnatural experiments. Causality is usually inferred by comparing the effect of a policy – or more generally of a treatment – with its counterfactual, ideally given by a random control group. The causal effect of the treatment is then defined as $\Delta=y_1-y_0$, i.e. the difference between what the situation would be with treatment (noted t=1) and without treatment (noted t=0). The concern is that only $y=t\cdot y_1+(1-t)\cdot y_0$ and $t$ are observed. In other words, the causal effect of variable $t$  on $t$  is not observed (since only one of the two potential variables – $y_0$ or $y_1$  is observed for each individual), but it is also individual, and therefore a function of x-covariates. Generally, by making assumptions about the distribution of the triplet $(Y_0,Y_1,T)$, some parameters of the causal effect distribution become identifiable, based on the density of the observable variables $(Y,T)$. Classically, we will be interested in the moments of this distribution, in particular the average effect of treatment in the population, $\mathbb{E}[\Delta]$, or even just the average effect of treatment in the case of treatment $\mathbb{E}[\Delta|T=1]$. If the result $(Y_0,Y_1)$ is independent of the processing access variable $T$, it can be shown that $\mathbb{E}[\Delta]=\mathbb{E}[Y|T=1]- \mathbb{E} [Y|T=0]$. But if this independence hypothesis is not verified, there is a selection bias, often associated with $\mathbb{E}[Y_0|T=1]- \mathbb{E} [Y_0|T=0]$. Rosenbaum & Rubin (1983) propose to use a propensity to be treated score, $p(x)=\mathbb{P}[T=1|X=x]$, noting that if variable $Y_0$\ is independent of access to treatment $T$ conditionally to the explanatory variables $X$, then it is independent of $T$  conditionally to the score $p(X)$ : it is sufficient to match them using their propensity score. Heckman et al (2003) thus proposes a kernel estimator on the propensity score, which simply provides an estimator of the effect of the treatment, provided that it is treated.

To be continued next time, we’ll introduce “machine learning techniques” (references mentioned above are online here)

# Probabilistic Foundations of Econometrics, part 3

This post is the third one of our series on the history and foundations of econometric and machine learning models. Part 2 is online here.

## Exponential family and linear models

The Gaussian linear model is a special case of a large family of linear models, obtained when the conditional distribution of $Y$ (given the covariates) belongs to the exponential family$$f(y_i|\theta_i,\phi)=\exp\left(\frac{y_i\theta_i-b(\theta_i)}{a(\phi)}+c(y_i,\phi)\right)$$ with $\theta_i=\psi(\mathbf{x}_i^T \beta)$. Functions $a$, $b$ and $c$ are specified according to the type of exponential law (studied extensively in statistics since Darmoix (1935), as Brown (1986) reminds us), and $\psi$ is a one-to-one mapping that the user must specify. Log-likelihood then has a simple expression $$\log\mathcal{L}(\mathbf{\theta},\phi|\mathbf{y}) =\frac{\sum_{i=1}^ny_i\theta_i-\sum_{i=1}^nb(\theta_i)}{a(\phi)}+\sum_{i=1}^n c(y_i,\phi)$$ and the first order condition is then written$$\frac{\partial \log \mathcal{L}(\mathbf{\theta},\phi|\mathbf{y})}{\partial \mathbf{\beta}} = \mathbf{X}^T\mathbf{W}^{-1}[\mathbf{y}-\widehat{\mathbf{y}}]=\mathbf{0}$$based on Müller’s (2011) notations, where $\mathbf{W}$ is a weight matrix (which depends on $\beta$). Given the link between $\theta$ and the expectation of $Y$, instead of specifying the function $\psi(\cdot)$, we will tend to specify the link function $g(\cdot)$ defined by $$\widehat{y}=m(\mathbf{x})=\mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=g^{-1} (\mathbf{x}^T \beta)$$For the Gaussian linear regression we consider an identity link, while for the Poisson regression, the natural link (called canonical) is the logarithmic link. Here, as $\mathbf{W}$ depends on $\beta$ (with $\mathbf{W}=diag(\nabla g(\widehat{\mathbf{y}})Var[\mathbf{y}])$ there is generally no explicit formula for the maximum likelihood estimator. But an iterative algorithm makes it possible to obtain a numerical approximation. By setting $$\mathbf{z}=g(\widehat{\mathbf{y}})+(\mathbf{y}-\widehat{\mathbf{y}})\cdot\nabla g(\widehat{\mathbf{y}})$$ corresponding to the error term of a Taylor development in order 1 of $g$, we obtain an algorithm of the form$$\widehat{\beta}_{k+1}=[\mathbf{X}^T \mathbf{W}_k^{-1} \mathbf{X}]^{-1} \mathbf{X}^T \mathbf{W}_k^{-1} \mathbf{z}_k$$By iterating, we will define $\widehat{\beta}=\widehat{\beta}_{\infty}$, and we can show that – with some additional technical assumptions (detailed in Müller (2011)) – this estimator is asymptotically Gaussian, with $$\sqrt{n}(\widehat{\beta} -\beta)\overset{\mathcal{L}}{\rightarrow} \mathcal{N}(\mathbf{0},I(β)^{-1})$$where numerically $I(\beta)=\varphi\cdot[\mathbf{X}^T \mathbf{W}_\infty^{-1} \mathbf{X}]$.

From a numerical point of view, the computer will solve the first-order condition, and actually, the law of $Y$ does not really intervene. For example, one can estimate a “Poisson regression” even when observations are not integers (but they need to be positive). In other words, the law of $Y$ is only an interpretation here, and the algorithm could be introduced in a different way (as we will see later on), without necessarily having an underlying probabilistic model.

## Logistic Regression

Logistic regression is the generalized linear model obtained with a Bernoulli’s law, and a link function which is the quantile function of a logistic law (which corresponds to the canonical link in the sense of the exponential family). Taking into account the form of Bernoulli’s law, econometrics proposes a model for $y_i\in\{0,1\}$, in which the logarithm of the odds follows a linear model: $$\log\left(\frac{\mathbb{P}[Y=1\vert \mathbf{X}=\mathbf{x}]}{\mathbb{P}[Y\neq 1\vert \mathbf{X}=\mathbf{x}]}\right)=\beta_0+\mathbf{x}^T\beta$$or $$\mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=\mathbb{P}[Y=1|\mathbf{X}=\mathbf{x}]=\frac{e^{\beta_0+\mathbf{x}^T\beta}}{1+ e^{\beta_0+\mathbf{x}^T\beta}}=H(\beta_0+\mathbf{x}^T\beta)$$where $H(\cdot)=\exp(\cdot)/(1+exp(\cdot))$ is the cumulative distribution function of the logistic law. The estimation of $(\beta_0,\beta)$ is performed by maximizing the likelihood: $$\mathcal{L}=\prod_{i=1}^n \left(\frac{e^{\mathbf{x}_i^T\mathbf{\beta}}}{1+e^{\boldsymbol{x}_i^T\mathbf{\beta}}}\right)^{y_i}\left(\frac{1}{1+e^{\mathbf{x}_i^T\mathbf{\beta}}}\right)^{1-y_i}$$ It is said to be a linear models because isoprobability curves here are the parallel hyperplanes $b+\mathbf{x}^T\beta$. Rather than this model, popularized by Berkson (1944), some will prefer the probit model (see Berkson, 1951), introduced by Bliss (1934). In this model: $$\mathbb{E}[Y|\mathbf{X}=\mathbf{x}]=\mathbb{P}[Y=1|\mathbf{X}=\mathbf{x}]=\Phi (\beta_0+\mathbf{x}^T\beta)$$

where $\Phi$ denotes the distribution function of the reduced centred normal distribution. This model has the advantage of having a direct link with the Gaussian linear model, since $y_i=\mathbf{1}(y_i^\star>0)$ with $y_i^\star=\beta_0+\mathbf{x}^T \beta+\varepsilon_i$ where the residuals are Gaussian, $\mathcal{N}(0,\sigma^2)$. An alternative is to have centered residuals of unit variance, and to consider a latent modeling of the form $y_i=\mathbf{1}(y_i^\star>\xi)$ (where $\xi$ will be fixed). As we can see, these techniques are fundamentally linked to an underlying stochastic model. In the body of the article, we present several alternative techniques – from the learning literature – for this classification problem (with two classes, here $0$ and $1$).

## Regression in high dimension

As we mentioned earlier, the first order condition $\mathbf{X}^T (\mathbf{X}\widehat{\beta}-\mathbf{y})=\mathbf{0}$ is solved numerically by performing a QR decomposition, at a cost which consists in $O(np^2)$ operations (where $p$ is the rank of $\mathbf{X}^T \mathbf{X}$). Numerically, this calculation can be long (either because $p$ is large or because $n$ is large), and a simpler strategy may be to sub-sample. Let $n_s\ll n$, and consider a sub-sample size $n_s$ of $\{1,\cdots,n\}$. Then $\widehat{\beta}_s=(\mathbf{X}_s^T \mathbf{X}_s )^{-1} \mathbf{X}_s^T\mathbf{y}_s$ is a good approximation of $\beta$ as shown by Dhillon et al. (2014). However, this algorithm is dangerous if some points have a high leverage (i.e. $L_i=\mathbf{x}_i(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{x}_i^T$). Tropp (2011) proposes to transform the data (in a linear way), but a more popular approach is to do non-uniform sub-sampling, with a probability related to the influence of observations (defined by $I_i=\widehat{\varepsilon}_iL_i/(1-L_i)^2$, and which unfortunately can only be calculated once the model is estimated).

In general, we will talk about massive data when the data table of size does not fit in the RAM memory of the computer. This situation is often encountered in statistical learning nowadays with very often $p\ll n$. This is why, in practice, many libraries of algorithms assimilated to machine learning use iterative methods to solve the first-order condition. When the parametric model to be calibrated is indeed convex and semi-differentiable, it is possible to use, for example, the stochastic gradient descent method as suggested by Bottou (2010). This last one allows to free oneself at each iteration from the calculation of the gradient on each observation of our learning base. Rather than making an average descent at each iteration, we start by drawing (without replacement) an observation $\mathbf{x}_i$ among the $n$ available. The model parameters are then corrected so that the prediction made from $\mathbf{x}_i$ is as close as possible to the true value $y_i$. The method is then repeated until all the data have been reviewed. In this algorithm there is therefore as much iteration as there are observations. Unlike the gradient descent algorithm (or Newton’s method) at each iteration, only one gradient vector is calculated (and no longer $n$). However, it is sometimes necessary to run this algorithm several times to increase the convergence of the model parameters. If the objective is, for example, to minimize a loss function $\ell$ between the estimator $m_\beta (\mathbf{x})$ and $y$ (like the quadratic loss function, as in the Gaussian linear regression) the algorithm can be summarized as follows:

• Step 0: Mix the data
• Iteration step: For $t=1,\cdots, n$, we pull $i\in\{1,\cdots,n\}$ without replacement, and we set $$\beta^{t+1} = \beta^{t} - \gamma_t\frac{ \partial{\ell(y_i,m_{\beta^t}(X_i)) } }{ \partial{ \beta}}$$

This algorithm can be repeated several times as a whole depending on the user’s needs. The advantage of this method is that at each iteration, it is not necessary to calculate the gradient on all observations (more sum). It is therefore suitable for large databases. This algorithm is based on a convergence in probability towards a neighborhood of the optimum (and not the optimum itself).

(references will be given in the very last post of that series) To be continued