{"word_count": 2193, "character_count": 20395}
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’sWe have a training sample, with observations [latex](\mathbf{x}_i,y_i)[/latex] where the variables [latex]y[/latex] are in a set [latex]\mathcal{Y}[/latex]. In the case of classification, [latex]\mathcal{Y}=\{-1,+1\}[/latex], but a relatively general set can be considered (note that if econometricians prefer [latex]\mathcal{Y}=\{0,1\}[/latex] - because of the Bernoulli distribution and because [latex]0[/latex] and [latex]1[/latex] are lower and upper bounds of probabilities, people in the "machine learning" community prefer [latex]\mathcal{Y}=\{-1,+1\}[/latex]). A predictor [latex]m[/latex] is an function taking values in [latex]\mathcal{Y}[/latex], used to label (or classify) future new observations, using some features that lie in a set [latex]\mathcal{X}[/latex]. It is assumed that the labels are produced by an (unknown) classifier [latex]f[/latex] called target. For a statistician, this function would be the real model. Naturally, we want to build [latex]m[/latex] as close as possible to [latex]f[/latex]. Let [latex]\mathbb{P}[/latex] be a (unknown) distribution on [latex]\mathcal{X}[/latex]. The error of [latex]m[/latex] with respect to target [latex]f[/latex] is defined by [latex display="true"]\mathcal{R}_{\mathbb{P},f}(m)=\mathbb{P}[m(\boldsymbol{X})\neq f(\boldsymbol{X})]\text{ where }\boldsymbol{X}\sim\mathbb{P}[/latex]or equivalently,[latex display="true"]\mathcal{R}_{\mathbb{P},f}(m)=\mathbb{P}\big[\{\boldsymbol{x}\in\mathcal{X}:m(\boldsymbol{x})\neq f(\boldsymbol{x})\}\big][/latex]To obtain our “optimal” classifier, it becomes necessary to assume that there is a link between the data in our sample and the pair [latex] (\mathbb{P},f) [/latex], i.e. a data generation model. We will then assume that the [latex]\mathbf{x}_i[/latex] are obtained by independent draws according to [latex]\mathbb{P}[/latex], and that then [latex]y_i=f(\mathbf{x}_i) [/latex] . We can define the empirical risk of a classifier [latex]m[/latex], as [latex display="true"]\widehat{{R}}(m)=\frac{1}{n}\sum_{i=1}^n \boldsymbol{1}(m(\boldsymbol{x}_i)\neq y_i)[/latex]It is important to recognize that a perfect model cannot be found, in the sense that [latex]R_{\mathbb{P},f} (m)=0[/latex]. Indeed, if we consider the simplest case, with [latex]\mathcal{X}=\{x_1,x_2\}[/latex] and [latex]\mathbb{P}[/latex] is such that [latex]\mathbb{P}(\{x_1\})=p[/latex] and [latex]\mathbb{P}(\{x_2\})=1-p[/latex]. The probability of never observing [latex]\{x_2\}[/latex] among the [latex]n[/latex] observations is [latex] (1-p)^n[/latex], and if [latex]p&#60;1/n[/latex], it is quite likely never to observe [latex]\{x_2\}[/latex] so it can never be predicted. We cannot therefore hope to have a zero risk whatever [latex]\mathbb{P}[/latex]. And more generally, it is also possible to observe [latex]\{x_1\}[/latex] and [latex]\{x_2\}[/latex], 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 [latex]m[/latex] such that [latex]R_{\mathbb{P},f} (m)\leq\varepsilon[/latex], where [latex]\varepsilon[/latex] is an a priori specified threshold. But even this condition is too strong, and cannot be fulfilled. Thus, we will usually as to have [latex]R_{\mathbb{P},f} (m)\leq\varepsilon[/latex] with some probability [latex]1-\delta[/latex]. Hence, we will try to be "probably approximately correct" (PAC), allowing to make a mistake with a probability [latex]\delta[/latex], again fixed a priori.Also, when we build a classifier, we do not know either [latex]\mathbb{P}[/latex] or [latex]f[/latex], but we give ourselves a precision criterion [latex]\varepsilon [/latex], and a confidence parameter [latex] \delta[/latex], and we have [latex]n[/latex] observations. Note that [latex]n[/latex], [latex] \varepsilon[/latex] and [latex]\delta[/latex] can be linked. We then look for a model [latex]m[/latex] such that [latex]R_{\mathbb{P},f} (m)\leq\varepsilon[/latex] with probability (at least) [latex]1-\delta[/latex], so that we are probably approximately correct. Wolpert (1996) has shown (see details in Wolpert &#38; Macready (1997)) that there is no universal learning algorithm. In particular, it can be shown that there is [latex]\mathbb{P}[/latex] such that [latex]R_{\mathbb{P},f} (m) [/latex] 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 [latex]m[/latex], we will then force [latex]m[/latex] to belong to a particular class, noted [latex]\mathcal{M}[/latex]. Let us suppose, to start with, that [latex]\mathcal{M}[/latex] contains a finite number of possible models. We can then show that for all [latex]\varepsilon[/latex] and [latex]\delta[/latex], that for all [latex]\mathbb{P}[/latex] and [latex]f[/latex], if we have enough observations (more precisely [latex]n\geq \varepsilon^{-1} \log[\delta^{-1} |\mathcal{M}|][/latex], then with a greater probability than [latex]1-\delta[/latex], [latex]R_{\mathbb{P},f} (m^\star)\leq\varepsilon[/latex] where[latex display="true"] 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[/latex]in other words [latex]m^\star[/latex] is a model in [latex]\mathcal{M}[/latex] that minimizes empirical risk.We can go a little further, staying in the case where [latex]\mathcal{Y}=\{-1,+1\}[/latex]. An [latex]\mathcal{M}[/latex] class of classifiers will be called PAC-learnable if there is [latex]n_M:[0,1]^2\rightarrow \mathbb{N}[/latex] such that, for all [latex]\varepsilon[/latex], [latex]\delta[/latex], [latex]\mathbb{P}[/latex] and if it is assumed that the target [latex]f[/latex] belongs to [latex]\mathcal{M}[/latex], then using [latex]n&#62;n_M (\varepsilon,\delta) [/latex] observations [latex]\mathbf{x}_i[/latex] drawn from [latex]\mathbb{P}[/latex], labelled [latex]y_i[/latex] by [latex]f[/latex], then there is [latex]m\in\mathcal{M}[/latex] such that, with probability [latex]1-\delta[/latex], [latex]R_{\mathbb{P},f} (m)\leq\varepsilon[/latex]. The [latex]n_M[/latex] function is then called "sample complexity to learn". In particular, we have seen that if [latex]M[/latex] contains a finite number of classifiers, then [latex]\mathcal{M}[/latex] is PAC-learnable with complexity [latex]n_M (\varepsilon,\delta)=\varepsilon^{-1} \log[\delta^{-1} |M|][/latex].Naturally, we would like to have a more general result, especially if [latex]\mathcal{M}[/latex] is not finite. To do this, the [latex]VC[/latex] dimension of Vapnik-Chervonenkis must be used, which is based on the idea of shattering points (for a binary classification). Consider [latex]k[/latex] points [latex]\{x_1,\cdot,x_k\}[/latex], and consider the set [latex display="true"]{E}_k=\big\lbrace(m(\boldsymbol{x}_1),\cdots,m(\boldsymbol{x}_k))\text{ for }m\in\mathcal{M})\big\rbrace[/latex] Note that the elements of [latex]E_k[/latex] belong to [latex]\{-1,+1\}^k[/latex]. In other words, [latex]|E_k |\leq 2^k[/latex]. We will say that M shatter all the points if all the combinations are possible, i. e. [latex]|E_k |=2^k[/latex]. Intuitively, the labels of the set of points do not provide enough information on target [latex]f[/latex], because anything is possible. The [latex]VC[/latex] dimension of [latex]\mathcal{M}[/latex] is then[latex display="true"] VC(\mathcal{M})=\sup\big\lbrace k\text{ such that }\mathcal{M}\text{ shatters }\{\boldsymbol{x}_1,\cdots\boldsymbol{x}_k\}\big\rbrace[/latex]For example, if [latex] \mathcal{X}=\mathbb{R}[/latex] and all (simple) models of the form [1] [latex]m_{a,b}=\mathbf{1}_{\pm}(x\in[a,b])[/latex] are considered. No set of [latex]\{x_1,x_2,x_2,x_3\}[/latex] ordered points can be shattered because it is sufficient to assign respectively +1, -1 and +1 to [latex]x_1[/latex], [latex]x_2[/latex] and [latex]x_3[/latex] respectively, therefore [latex]VC&#60;3[/latex]. On the other hand [latex]\{0,1\}[/latex] is shattered, so [latex]VC\geq 2[/latex]. The dimension of this predictor set is [latex]2[/latex]: If we increase by one dimension, [latex]\mathcal{X}=\mathbb{R}^2[/latex] and consider all (simple) models of the form [latex]m_{a,b}=\mathbf{1}_{\pm} (x\in[a,b]) [/latex] (where [latex] [a,b] [/latex] refers to the rectangle), then the dimension of [latex]\mathcal{M}[/latex] is here [latex]4[/latex].To introduce SVMs, let's place ourselves in the case where [latex]\mathcal{X}=\mathbb{R}^k[/latex], and consider separations by hyperplanes passing through the origin (we will say homogeneous), in the sense that [latex]m_{\mathbf{w}} (\mathbf{x})=\mathbf{1}_{\pm}(\mathbf{w}^T \mathbf{x}\geq 0) [/latex]. It can be shown that no set of [latex]k+1[/latex] points can be shattered by these two homogeneous spaces in [latex]\mathbb{R}^k[/latex], and therefore [latex]VC(M)=k[/latex]. If we add a constant, in the sense that [latex]m_{\mathbf{w},b} (\mathbf{x})=\mathbf{1}_{\pm}(\mathbf{w}^T \mathbf{x}+b\geq 0)[/latex], we can show that no set of [latex]k+2[/latex] points can be sprayed by these two (non-homogeneous) spaces in [latex]\mathbb{R}^k[/latex], and therefore [latex]VC(M)=k+1[/latex]. This dimension reminds us of the dimension of the model we've seen in the econometric context.From this dimension [latex]VC[/latex], we deduce the so-called fundamental theorem of learning: if [latex]\mathcal{M}[/latex] is a class of dimension [latex]d=VC(M) [/latex], then there are positive constants [latex]\underline{C}[/latex] and [latex]\overline{C}[/latex] such as the sample complexity for M to be PAC-learnable satisfies[latex display="true"] \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)[/latex]The link between the notion of learning (as defined in Vailiant (1984)) and the [latex]VC[/latex] 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 &#38; Hart (1967) on learning in the context of the algorithm of the [latex]k[/latex]-nearest neighbors. These studies have linked learning, information theory (with the textbook Cover &#38; 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 functionThese 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 [latex]y\in\{-1,+1\}[/latex], inspired by McCulloch &#38; Pitts (1943). We have data [latex]\{(y_i,\mathbf{x}_i)\}[/latex], and we will iteratively build a set of [latex]m_k[\mathbf{x}[/latex] 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 :[latex display="true"]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}&#60; 0\end{array}\right.[/latex]where [latex]\beta[/latex] coefficients are often interpreted as "weights" assigned to each of the explanatory variables. We give ourselves initial weights [latex] (\beta_0^{(0)},\beta^{(0)} [/latex], which we will update taking into account the prediction error made, between [latex]y_i[/latex] and the prediction [latex]\widehat{y}_i^{(k)}[/latex] :[latex display="true"]\widehat{y}_i^{(k)}=m^{(k)}(\mathbf{x}_i)=\boldsymbol{1}_{\pm}(\beta_0^{(k)}+\mathbf{x}^T \boldsymbol{\beta}^{(k)}\geq 0), [/latex]with, in the case of the perceptron:[latex display="true"]\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[/latex]Here [latex]\ell(y,y')=\mathbf{1}(y\neq y') [/latex] is a loss function, which will allow to give a price to an error made, by predicting [latex]\widehat{y}=m(\mathbf{x})[/latex] and observing [latex]y[/latex]. For a regression problem, we can consider a quadratic error [latex]\ell_2[/latex], such that [latex] \ell(y,m(\mathbf{x}))=(y-m(\mathbf{x}))^2[/latex] or in absolute value [latex]\ell_1[/latex], with [latex] \ell(y,m(\mathbf{x}))=|y-m(\mathbf{x})|[/latex]. 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:[latex display="true"]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)[/latex]for a predefined set of predictors [latex]\mathcal{M}[/latex]. Any machine learning problem is mathematically formulated as an optimization problem, whose solution determines a set of model parameters (if the [latex]\mathcal{M}[/latex] family is described by a set of parameters - which can be coordinates in a functional database). We can note [latex]\mathcal{M}_0[/latex] the space of the hyperplanes of [latex]\mathbb{R}^p[/latex] in the sense that[latex display="true"]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[/latex] 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 [latex]\widehat{m}^\star[/latex], known as "oracle", in a family of [latex]\mathcal{M}[/latex] estimators, [latex display="true"]\widehat{m}^{\star} =\underset{\widehat{m}\in\mathcal{M}}{\text{argmin}}\big\lbrace\mathcal{R}(\widehat{m},m)\big\rbrace[/latex]This estimator is, of course, impossible to define because it depends on [latex]m[/latex], the real model, unknown.But let's come back a little more to these loss functions. A loss function [latex]\ell[/latex] is a function [latex]\mathbb{R}^d\times\mathbb{R}^d\rightarrow\mathbb{R}_+[/latex], symmetric, which checks the triangular inequality, and such that [latex]\ell(x,y)=0[/latex] if and only if [latex]x=y[/latex]. The associated norm is [latex]\|\cdot\|[/latex], such that [latex]\ell(x,y)=\|x-y\|=\ell(x-y,0)[/latex] (using the fact that [latex]\ell(x,y+z)=\ell(x-y,z)[/latex] - 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:[latex display="true"]\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[/latex] where [latex]\ell_2[/latex] is the usual quadratic distance If we assume - as we did in econometrics - that there is an underlying probabilistic model, and observe that : [latex display="true"]\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}[/latex]it should be noted that what we are trying to obtain here, by solving the problem [latex](6)[/latex] by taking the norm [latex]\ell_2[/latex], is an approximation (in a given functional space, [latex]\mathcal{M}[/latex]) of the conditional expectation [latex]x\mapsto\mathbb{E}[Y|\mathbf{X}=\mathbf{x}][/latex]. Another particularly interesting loss function is the loss [latex]\ell_1[/latex],[latex] \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}[/latex]The optimization problem :[latex display="true"]\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[/latex] is obtained in econometrics by assuming that the conditional law of [latex]Y[/latex] follows a Laplace law centered on [latex]m(\mathbf{x})[/latex], 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 [latex]Y[/latex] is symmetrical with respect to [latex]0[/latex], the median and the mean coincide If this loss function is rewritten  [latex display="true"] \ell_1(y,m)=\vert (y-m)(1/2-\boldsymbol{1}_{y\leq m})\vert[/latex] a generalization can be obtained for [latex]\tau\in[0.1][/latex]:[latex display="true"]\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[/latex]where[latex display="true"]\ell_{\tau}^{q}(x,y)= (x-y)(\tau-\boldsymbol{1}_{x\leq y})[/latex]  is then the quantile regression of level [latex]\tau[/latex] (Koenker, 2003; d'Haultefœuille &#38; 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: [latex display="true"]\displaystyle{\ell}^{\text{ e}}_{\tau}(x,y)= (x-y)^2\cdot\big\vert\tau-\boldsymbol{1}_{x\leq y}\big\vert[/latex]with [latex]\tau\in[0.1][/latex]. We see the parallel with the quantile function: [latex display="true"]\displaystyle{\ell}^{\text{ q}}_{\tau}(x,y)= \vert x-y\vert \cdot\big\vert\tau-\boldsymbol{1}_{x\leq y}\big\vert[/latex]Koenker &#38; Machado (1999) and Yu &#38; Moyeed (2001) also noted a link between this condition and the search for maximum likelihood when [latex]Y[/latex]'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 [latex]T[/latex] will be said to be "ellicitable" if there is a loss function [latex]\ell:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}_+[/latex] such that:[latex display="true"]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[/latex] The mean (mathematical expectation) is thus ellicable by the quadratic distance, [latex]\ell_2[/latex], while the median is ellicable by the distance [latex]\ell_1[/latex]. 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 &#38; 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 [latex]\mathbf{1}_{\pm}[/latex] does not take values 0 or 1 (like the classical [latex]\mathbf{1}[/latex] function), but -1 and +1.