Machine Learning 1. Lecture 12. Density Estimation and Bayesian Classification

What this lecture is about

The course programme actually has two lectures here: one on density estimation on its own, the other on Bayesian classification theory on its own. They are closely related things, so this semester I merged them, but the merged lecture bears the traces of the merger: it has come out large, and some things I will tell in less detail or leave out. That should not affect the general understanding.

Today it is the problem of estimating a probability density, and we will look at three main approaches: the classical parametric one, mixtures of distributions, and the nonparametric one. Then, on the basis of these approaches, we will see how classification algorithms are built.

Parametric density estimation

The maximum likelihood principle

Suppose we have a sample. This is unsupervised learning: there are objects, but there are no answers for them, only $x$. It is assumed that these $x$ came from some density $p(x)$, and this density has to be recovered: we have to understand which density the observed objects of the sample came from.

The first thing that comes to mind: let us describe the density by a parametric family of functions. Say we know the form of the dependence, but the parameters are unknown. This is very much like the way we solve classification and regression problems: we fix a dependence with parameters and train the parameters. But we have to understand by what criterion to train, and here mathematical statistics comes to our aid, telling us that such a problem should be solved by the maximum likelihood principle. We write out the probability density of the whole sample: it is the product of the one-object densities. We are entitled to write a product because there is the i.i.d. assumption, independent identically distributed, a simple sample: the observations appear independently of one another from one and the same density. The logarithm is a monotone function, so we may maximise the logarithm. And if we are used to minimising a loss function, then our loss function is minus the logarithm of the density model. The problem has again been reduced to optimisation, and numerical methods can be used. But let me tell you about the classical situation where this problem can be solved analytically.

The multivariate Gaussian density

This is the multivariate Gaussian distribution. Why the density formula is like this is material for another course, mathematical statistics; for us it is a working tool. The function is given up to its parameters: $\mu$ is the expectation vector, $\Sigma$ the covariance matrix. We have to learn to estimate them from the sample, of course by the maximum likelihood principle. For that the function has to be carefully differentiated with respect to the parameters, $\mu$ and $\Sigma$. $\mu$ is a vector, and we know how to differentiate with respect to a vector, but $\Sigma$ is a covariance matrix, a symmetric one at that, so differentiating with respect to it is not a very simple undertaking. There is an inverse matrix in there, and the determinant of $\Sigma$ in the denominator, and one has to learn to differentiate a scalar function, the density, with respect to such expressions. We will skip these technical details now; the proofs are textbook ones and are given in many places. The main thing is that the result is very clear and predictable.

The sample estimate of the expectation vector is, of course, the mean over all the objects of the sample: the centre of the sample is easy to compute, we just average all the vectors $x_i$. And the sample estimate of the covariance matrix is, as one would expect, this matrix: $x_i-\hat{\mu }$ are the centred object vectors, and from each of them we made a matrix by multiplying a column by a row, the vector by itself. The sum of $\ell$ such matrices is the sample estimate of the covariance matrix.

The geometric meaning

What a multivariate Gaussian density looks like you should also already understand from other courses; this is just a reminder. It is a swarm of points, an ellipse; it is called the scatter ellipsoid. One can show that the axes of this ellipse are connected with the orthogonal eigenvectors of the covariance matrix. We will not do that now; today's lecture does not use it specifically. But it is useful to illustrate the geometric meaning of the multivariate Gaussian, simply to picture it as a cloud of points of elliptical shape in the space ${\mathbb{R}}^n$, the space of object vectors.

The multicollinearity problem

There is a little problem, of course: we are going to invert the covariance matrix, and it may turn out to be singular, and in some cases it will be singular with a guarantee. Let us look again at how the sample estimate of the covariance matrix is defined: it is the sum of $\ell$ matrices of rank one. Imagine the situation where the amount of data $\ell$ is less than the dimension of the vectors, that is, the number of features $n$. You have added up $\ell$ matrices of rank one, and the rank of the resulting matrix cannot exceed $\ell$; it is of size $n\times n$, and there are not enough objects. Then the estimate of the covariance matrix is singular with a guarantee; it is a matrix of incomplete rank. And even if there are more objects than features, it may still turn out to be close to singular: there may be linear dependences among the features, multicollinearity. We have already discussed how to work with such things when we talked about regression. A slight modification of the covariance matrix is possible that at once makes it well-conditioned; one only has to choose the constant $\tau$, the regularisation coefficient.

Let us look at the scatter ellipsoid again. Adding $\tau$ to the diagonal leads to every eigenvalue increasing by exactly $\tau$, while the eigenvectors stay where they are: the scatter ellipsoid gets inflated a little, along every coordinate it grows by $\tau$. What is multicollinearity in this geometric interpretation? It is an ellipse flattened along some coordinates right down to zero. And when you have inflated it by adding $\tau$ to each eigenvalue, the covariance matrix has ceased to be singular. That is the geometric meaning of regularisation: inflating the scatter ellipsoid.

Or one can do something quite radical and say: let the covariance matrix be diagonal altogether. On its diagonal are the variances of the features, and off the diagonal the covariances between features. If we assume that the features are independent, and for Gaussian distributions this is equivalent to their being uncorrelated, the problem becomes very much simpler: we no longer need to invert non-diagonal covariance matrices, and diagonal ones are easily inverted. Estimating an $n$-dimensional density becomes a very simple problem: estimating $n$ one-dimensional densities. For each feature separately we compute the expectation and the variance: the ordinary familiar scalar estimates, the mean value of the feature over all objects and its variance. Of course, such an assumption rarely holds: in most problems one should hardly expect the features to be statistically independent. But the solution becomes much simpler: recovering $n$ one-dimensional densities is much easier than recovering one $n$-dimensional density.

That is a short story about how density estimation problems are solved if we know the functional form of the density. Analytically the problem is solved neatly only for the normal density; for more complex forms of densities one has to use numerical optimisation.

Separating a mixture of distributions

Now let us consider another situation where there is a numerical method, but a specific one, very interesting and used in many places in machine learning. We will look at it now for separating a mixture of distributions. You have some data sample and do not believe that it is represented by a single scatter ellipsoid; you cannot describe it by one multivariate Gaussian density. But you suspect that there is a mixture of densities there: the complex shape of the sample is due to part of the data coming from one Gaussian density, part from another, part from a third. Though why Gaussian? Let us do it in general form.

A mixture of distributions

There are separate components $\phi$, which are still probability densities, but there are several of them, and their probabilistic mixture with coefficients that are in essence a discrete probability distribution. The expression for a mixture of distributions describes a two-level generative model of the data. What does two-level mean? $x$ is generated in two steps. First, from the discrete distribution of the prior probabilities of the components, a kind of asymmetric coin or $k$-sided die, you toss it and generate the number $j$ of a mixture component. Then you go to this component and generate $x$ itself from its distribution. Such two-stage schemes for generating objects are also called a generative story: a little tale of how the objects of the sample appear. They no longer come from one distribution but from a mixture of distributions.

A mixture of distributions is in fact also a parametric form of writing a distribution, only it has more parameters: now we have to estimate the parameters ${\theta }_j$ of each mixture component, and there are also the non-negative normalised prior probabilities of the components $w_j$; they have to be estimated too. If we write out the problem of maximising the log-likelihood, we get a very inconvenient functional, where there is a sum under the logarithm. Why was the problem solved analytically for the multivariate Gaussian distribution? Because it is exponential: the logarithm and the exponential immediately cancelled each other, and what was left was an expression that is easy to differentiate with respect to the model parameters, to get a system of equations and solve it. Here the sum under the logarithm at once cuts off every path to an analytical solution. Whichever way you turn it, a numerical method will have to be invented. One could, of course, solve this problem head-on with some stochastic gradient, but mankind has come up with a much more interesting way.

The EM algorithm

This way is called the EM algorithm, from the words expectation and maximization; that is what the two steps of this algorithm were called. It is an iterative process, and we will even prove a theorem now.

The theorem gives necessary conditions for an extremum: the totality of all the parameters of the mixture at a point of local extremum of our functional, the log-likelihood, satisfies a system of equations in all the parameters. And the system of equations is written in such a fortunate way that it simply asks to be solved by the method of simple iteration: just compute the expressions on the left-hand side of the system by these formulas, variable by variable. We only have to explain what these auxiliary variables $g_{ij}$ are. Looking at the formula, we see only that it is some function of the model parameters given by a not very complicated formula; it can be computed if $w_j$ and ${\theta }_j$ are known. That is, we have some initial approximation, from it we compute $g_{ij}$, and then, knowing $g_{ij}$, we can compute an approximation for $\theta$ and for $w$. That is the whole method of solving this system of equations.

Let us interpret this system, and then deal with the proof. It turns out that the E-step is simply Bayes' formula, and the $g_{ij}$ have a very clear probabilistic interpretation: they are the posterior probabilities of the components for the objects of the training sample. The meaning is this: if we look at the $i$-th object, with what probability could it have been generated by each, the $j$-th, component? That is what we estimate, and those are our auxiliary variables. Why Bayes' formula? We simply write it out by the formula for conditional probability. Note that I denote probabilities by a capital $P$ and probability densities in the space $X$ by a small $p$; I always follow this mnemonic rule. The numerator turns out to be one term of our mixture of distributions and the denominator the mixture itself, and that is where the formula came from. So far this is only an interpretation; we will now prove this formula rigorously, and it will appear from the solution of the optimisation problem as if by itself.

And the M-step is the solution of a problem very similar to maximising the log-likelihood, written out separately for each component in terms of the hidden variables $g_{ij}$. These auxiliary variables show to what degree each object of the sample belongs to the $j$-th component. Look: in the log-likelihood functional we took each $i$-th object with a weight coefficient, the degree to which it belongs to the $j$-th component. These weights are normalised over $j$. If an object does not belong to the $j$-th component at all, it will be removed from this sum, and in the sum we take only those objects that, as we estimated, were generated by the $j$-th component. It turns out that the problem of estimating the parameters of a mixture of distributions has split into $k$ independent subproblems of estimating each component separately. The E-step estimated the degree to which an object belongs to each component, and then we know which objects each component consists of, and we estimate each component by maximum likelihood. That is, we have reduced the problem of estimating a mixture to a problem we can already solve after a fashion, and for Gaussian components can solve well.

Now the proof, literally in a couple of words, because we have no time today, and I want to tell you everything.

The theorem does not in fact speak of an algorithm: it says that there is a system of equations, the necessary conditions for an extremum of the functional being optimised. It just happened to be written so, we were so lucky, that it can be solved by the method of simple iteration, and that method is what is called the EM algorithm. There are deeper proofs that this algorithm converges, unfortunately in a weak sense; there are no strong results about convergence for the EM algorithm, still less about global convergence. But the approach has proved very productive.

Let us gather everything into one piece of pseudocode. There is some initialisation: $\theta$ we will probably initialise randomly, and for the prior probabilities of the components we take the uniform distribution. Then we repeat the E-step and the M-step until the parameters converge. At the E-step we estimate by Bayes' formula to what degree each $i$-th object belongs to the $j$-th component. The M-step is a decomposition: each component of the mixture is estimated separately.

A mixture of Gaussians (GMM)

If we substitute multivariate Gaussian distributions for $\phi$, we get a very well-known variety of the EM algorithm: separating a mixture of Gaussians, the Gaussian Mixture Model, GMM. In general form $\theta$ is estimated as the solution of the problem of maximising the weighted log-likelihood, and it can be solved only when the concrete density $\phi$ is known. In the Gaussian case, fortunately, this problem is solved analytically: here are the analytical estimates for the expectation vector and for the covariance matrix. The summation here is over the whole sample, but with the weights $g_{ij}$: the $j$-th mixture component consists, as it were, of all the objects, but each object with its own weight. Some objects have a weight near zero, which means they practically do not belong to this component and will practically not be taken into account, and some may have a large weight: those are the objects from this component.

Let us look at a picture of how this process goes geometrically; it always becomes clearer that way. On the slide there are two Gaussian components in ${\mathbb{R}}^2$, and their initial approximations were deliberately chosen none too well: a blue circle and a red circle. The sample is a mixture of two distributions; note that if I had tried to stretch one Gaussian over it, that would not have been very good. At the first iteration the E-step estimates by Bayes' formula the degree of membership in each of the components. The components are so far poorly defined, but the E-step has already said: here are the blue objects that are closer to this component, the value of the blue density on them is much larger than of the red, so they were assigned to the blue component. Next step: from these "more blue than red" objects I estimate the density, and this is where it has moved; from the "more red than blue" ones, the second component, and this is where it has moved. Already at the first iteration these two components, first, covered the data, and second, began to move apart, each towards its own cluster. After the second iteration they shifted more, after the fifth more still, and after the twentieth it seems they already cover their clusters perfectly.

A very powerful approach. What is the meaning of the EM algorithm? It allowed us to crack the internal structure of the data. We had points, we did not know how they were structured, but we made the assumption that this structure is a mixture of two Gaussians, and the EM algorithm showed where these Gaussians are actually located.

Nonparametric density estimation

Now the next piece: the problem of nonparametric density estimation. What is wrong with parametric estimation? That we have to know an adequate model of the density. Where do we know it from? It will not always be possible to describe everything in the world by Gaussians, their mixtures or some other parametric distributions. Sometimes one can see that the data are non-Gaussian, but how they are distributed is unknown. What to do? Nonparametric approaches were invented. They arose and were actively developed somewhere in the middle of the 20th century: parametric approaches had long been known, they are the classics, and when people ran into really difficult problems, they started thinking how to do without a parametric model of the density.

The one-dimensional case

Let us begin with what a density is at all, and with the simplest case, where the objects are one-dimensional and, moreover, discrete. How do we estimate the density in this case? By a histogram: we simply compute the frequency estimate of $x$ taking a given value, and get frequency estimates of the probabilities of a discrete distribution.

And in the continuous case? We recall the definition of density. If there is a probability measure with which segments can be measured, let us take a segment centred at the point $x$ of width $2h$, $h$ to the left and $h$ to the right, and divide the probability of falling into this segment by its length $2h$. The limit as the length of the segment tends to zero is the density. The definition is non-constructive: there is a limit here, and one has to estimate the density somehow through the probability. What can be done here? Since there is a probability here, let us make a frequency estimate of it: in the neighbourhood of each point $x$ where we want to estimate the density, we look at what share of the objects of the sample falls into this neighbourhood. This is a frequency estimate of the probability: instead of $P$ we have substituted the share of the sample's objects, and the condition in brackets is the condition that $x_i$ has fallen into the $h$-neighbourhood of the point $x$. And we must not forget to divide by $2h$. We took the mathematical but non-computable definition, substituted a frequency estimate for the probability, and removed the limit. Clearly, if $h$ is made to tend to zero, the estimate degenerates: there will be delta functions at the points of the sample and zeros outside them. So the window width can no longer be made to tend to zero: it is a parameter, and it has to be determined somehow.

The Parzen–Rosenblatt estimate

Next we have to learn to estimate continuous densities. The term we are summing is a piecewise-constant function: the condition that $x_i$ lies in the $h$-neighbourhood of the point $x$, as a function of $x$, is a step, one in the neighbourhood of $x_i$ and zero outside it. We sum such steps, and it is somehow not right to recover continuous densities in the class of piecewise-constant functions. Can we replace this step by an approximation of it, preferably a smooth one? We can. Let us introduce a function that will be called a kernel, and require of it that it be even and normalised, the integral over the whole space equal to one; then the estimate, integrated over $x$, will also give one. And that it not increase: the farther a point is from $x$, the smaller its weight in the density. That is logical. Three requirements, and what exactly the kernel should look like we will still think about.

The multivariate case

We keep generalising: we began with the histogram, then recalled the definition of a continuous density, then said that the density estimate will be continuous too. And now let us remember that in our problems $x$ are usually multidimensional. What to do?

The first option: let us describe the density in the neighbourhood of the point $x$ by a product of one-dimensional densities. That can be done, though it is not very good: it means that in the neighbourhood of this point we consider the features locally independent of one another. True, afterwards we add up all such local density estimates, and after that the assumption of independence of the features is no longer there. The approach is good, but it is complicated by the fact that a window width of its own will have to be chosen for each feature, and that is not very convenient.

Another view. What did we have in the one-dimensional case? The difference $x-x_i$, and the kernel function is even, so it makes no difference whether we put the difference or the absolute value of the difference there. And the absolute value of the difference is already the distance between the two points $x$ and $x_i$. Let us imagine that in the multivariate case we have some metric in the space $X$, and substitute it into the density formula. But in this way we have created a problem for ourselves: we have to ensure that the function is integrable and its integral over the whole space equals one, so that it is a density. So a normalising factor will appear here, a function of the window width $h$, and it will have to be estimated somehow, which may not be very convenient.

The good news: these two approaches, at first sight completely different, simply coincide in some special cases. One such special case is the spherical Gaussian kernel, a special case of both. The same window width $h$ is introduced for all features, which in effect plays the role of the variance in a Gaussian distribution, and there are no dependences between the features: locally, in the neighbourhood of each point, we consider the features independent. The products of one-dimensional Gaussians over all features are then averaged over all objects of the sample. By way of generalisations we have come from a very simple one-dimensional case to what is called the Parzen density estimate, or the Parzen–Rosenblatt estimate.

The choice of kernel

Which kernel to choose? The piecewise-constant, rectangular kernel is an unfortunate choice, because the density estimate will come out piecewise-constant. And all the others are approximations, including smooth ones, such as the Gaussian. What is wrong with the Gaussian? It takes into account all the objects of the sample; it is not a localised density. And what if we take, for example, the quartic kernel, a piece of a fourth-degree polynomial on the segment $[-1,1]$? What is good about it: first, the function is smooth, even at zero, unlike the piece of a parabola, which is continuous but not smooth at the ends of the segment; and second, it is finite: only on the segment $[-1,1]$ is it non-zero, and outside it it equals zero. The quartic kernel is a very frequent choice; it really is often used.

There are theoretical results showing that in fact the choice of kernel has almost no effect on the quality of density estimation, but it does, of course, affect the smoothness of the resulting density. I will not dwell long on this result: practically all kernels give a value of the quality criterion very close to the optimal one, and the optimal one, strangely enough, is given by a piece of a parabola that is called the Epanechnikov kernel. On the other hand, the criterion raises questions: why do we compare densities in effect by the method of least squares? That is not a very good thing to do, but as for why, let us leave it for now; I will tell you some day.

The choice of window width

But the width of the smoothing window is precisely the parameter on which the quality of density estimation depends critically. The smoothness of the density estimate depends on the kernel, and its quality on the window width. On the slide there is an example: the sample is one-dimensional, blue points on the horizontal axis, clearly Gaussian; in grey the graph shows the true density, a Gaussian with zero expectation and unit variance. We try to recover it with different window widths. With the optimal window width we get a more or less decent approximation: the black curve is similar to the grey one. Take the window too small, and the estimate starts jerking after every point of the sample; these little spikes, these outliers arise because the density comes out too heavily jagged. This is overfitting: we start fitting to particular data points. And if we take $h$ too large, there will be underfitting; as $h$ tends to plus infinity, the estimate is smoothed out and tends to a uniform distribution. There is always a golden mean between the extreme cases: a window that is too small gives a noisy, overfitted model; one that is too large gives a model that does not describe the features of the data.

The typical form of the dependence $\mathrm{LOO}(h)$ (slide): for small $h$, overfitting; for large ones, underfitting; the optimal window width is in the neighbourhood of the minimum.

This is probably the main understanding: in Parzen density estimates it is very important to choose the window width; the quality of density estimation depends critically on it, and not much on the form of the kernel. The window width is very often chosen by leave-one-out: we train on the whole sample except one single point, and at that point we look at the loss function, the log-likelihood with a minus sign. Remember: the maximum likelihood principle maximises, so if we minimise, we must put a minus in front of the logarithm. We get a dependence $\mathrm{LOO}(h)$ with a minimum: on the left, for small $h$, overfitting; on the right, for large ones, underfitting. We look for the neighbourhood of the minimum and say: here we have determined the optimal window width. That is what is done most often.

A retrospective: three problems with one construction

This is very similar to what we have already seen when we considered the Parzen window method for metric classification algorithms and the method of nonparametric regression, the Nadaraya–Watson formula. Let us look at them all together and realise what we have just done. What is a nonparametric density estimate? It is a metric method. "Nonparametric" and "metric" are a play on words: they seem to be different things. "Nonparametric" means that we do not use a parametric model, and "metric" that we use a distance function between objects. At first sight these are quite unrelated things, but here they almost coincide, and it is correct to use either name.

We looked at the Parzen window method as a heuristic generalising the nearest neighbour method, and at its basis lay the distance between objects. Then we saw nonparametric regression methods, the Nadaraya–Watson kernel smoothing formula, and there too everything is based on the distance function between objects. And now we have realised that density estimation is the same kind of thing. All three problems have one common construction: a kernel and the choice of the width of the smoothing window in this kernel. The choice of kernel determines the smoothness of the function, whatever it is, a regression dependence, a separating surface in a classification problem or a density function, and the accuracy of the approximation is determined by the window width. This undoubtedly unites all three problems. The problems are different: classification and regression are supervised learning, density is unsupervised learning, but the approaches are very similar.

Three approaches to density estimation

Between the parametric approach, the nonparametric one and the separation of a mixture of distributions there is also an interesting fundamental connection. Look: nonparametric density estimation is in effect a mixture of distributions, with a distribution of its own at each data point. What is a kernel? It is a density in itself, and we have added up $\ell$ such densities. So this is a mixture, and with equal weights at that: each point, each object, has generated its own mixture component of some not very large width, and together in the mixture they have formed the density estimate.

You asked how to determine the optimal number of mixture components. And here I am showing that they can even be set equal to the number of objects of the training sample. True, such a thing cannot be tuned by the EM algorithm: there would be very heavy overfitting; a sample of $\ell$ objects cannot be described by $\ell$ mixture components. But here it turned out that it can: for each data point we built a mixture component of its own, with this point at its centre. This is a very extreme special case of a mixture, where the number of components is chosen not at all optimally, but since we fixed the window width, it does not lead to bad consequences such as overfitting, and the window width can be chosen for all mixture components at once. They are all of the same shape and width, and each is placed at its own object of the sample.

It turns out that recovering a mixture of distributions is, one might say, the golden mean between the ordinary parametric approach, which is a mixture of one component, and the nonparametric one, which is a mixture of $\ell$ components, where $\ell$ is the number of objects in the sample. And a mixture of $k$ components, if $k$ is determined optimally, is the most balanced approach: it is more general than both the parametric and the nonparametric one. That is an interpretation of how these three approaches to estimating probability densities are positioned: the mixture is probably the most general thing there can be.

Bayesian classification theory

Now let us take up classification; this is the second half of the lecture, though literally twenty minutes or so are left for it. Look what a very simple idea: if, when solving a classification problem, we learned to recover the density of each class, could we, knowing the estimates of the class densities, easily build a classifier? The answer is yes. This is called Bayesian classification theory: reducing the classification problem to the problem of recovering the distribution in each class separately.

Generative models of classification

Let us begin with the formulation. Now the objects are distributed over classes, and we consider that the joint space of "object–answer" pairs is a probability space with an unknown joint density of objects and answers, and our simple, i.i.d. sample came from this density. We are required to build an optimal classifier. Let us suppose for now that the joint density is known; later, of course, we will estimate it from the sample. If it were known, we would write it out by the formula for conditional density, or conditional probability, in two ways.

Let us look carefully at these expressions; here it is important not to get confused and to remember once and for all what is what, because this approach is used in many other places in machine learning. Of the four quantities, these are the ones we will need. The first is the prior probability of a class: what share of the objects of the sample lies in each class at all. Prior probabilities we can estimate by frequencies, the share of the sample in the class $y$. Next is the likelihood function of a class. What kind of beast is that? If we look at the class $y$, it is the density that the objects of this class have. This is exactly what we can recover as a density from the training sample: we take only the objects of the class $y$ and recover the density from them, which we have just learned to do. And the last, $P(y\mid x)$, is the posterior probability of a class: a posteriori because it appears after we have been given a fixed object $x$ and asked which class this $x$ could belong to. By Bayes' formula it is this thing.

Excellent. To which class, then, should the object $x$ be assigned? To the one whose posterior probability is maximal. The $\arg \max$ construction appears. But note: the $\arg \max$ is taken over the classes $y$, and the denominator $p(x)$ does not depend on $y$. It does not even need to be estimated; it can simply be thrown out: under the $\arg \max$ over $y$ it is just a constant. The joint density of all $x$ in all classes can be computed, but we simply do not need it. We need to estimate the prior probabilities of the classes and the densities of the objects in each of the classes, and this gives the Bayes classifier, so called because it is built by Bayes' formula.

An illustration: if the prior probabilities $P(y)$ were equal, the classification rule would become quite simple. At each object $x$ that we want to classify, we have to compute the density estimates for each of the classes and assign $x$ to the class for which the density estimate is maximal. Two classes, each described by its own density, and they have intersected; it seems that the intersection is the boundary between the classes, and we will need this consideration again. In one region one class has the larger density, in the other the other. A very simple visual interpretation of the Bayes classifier.

The discriminative and generative approaches

This is in fact a very deep consideration about how to solve classification problems in general: whether to build a separating surface or to recover the probability densities of the classes. It is a very serious dichotomy in classification methods, and the approaches are called discriminative and generative. The discriminative one says: at each object $x$ we will recover only the posterior probability, estimate with what probability the object belongs to a class. For this one can use, for example, a geometric approach: we build some separating surface, and the objects that turn out to be far from it belong to their class with high probability. This is like the support vector machine, which really is based on geometric assumptions; or logistic regression, where the posterior probability of the classes is modelled explicitly through a linear separating surface. The other approach is generative: we want to understand which densities the objects came to us from, to build density models for each class. This is in effect a model of data generation: each class had its own generating density, and we want to find out which. There are quite a lot of such examples too: for example, Parzen Window, the Parzen window method, or radial basis functions; the other abbreviations we do not know yet.

Which of the approaches is better cannot be said: sometimes one, sometimes the other, depending on the situation. But in the generative approach we are clearly doing some extra work: we try to describe the densities of the classes, including the objects that lie far from the separating surface. Why, one asks? We need only the classification, only the separating surface, which means the densities need to be described at the boundary between the classes, while we try to describe them entirely, over the whole space. It seems that we are doing extra work and hunting for peripheral objects that do not affect the quality of classification. That is why generative models need more data to train. But at the same time, since we recover the density, then, if the density model is decent, it will allow outliers to be ignored better: the robustness to noise is greater. When we discussed SVM, we said that the support objects turn out to be not only the objects on the boundary but also the objects that have jumped over the boundary, have become outliers, lie in another class, and we still use them as support objects. That is probably not very good. A generative model would tell us that these are objects at which the value of the density is very small, close to zero, so they can and should be ignored. There are many attempts to combine the merits of the discriminative and generative theories of classification.

The optimal Bayes classifier

I do not prove the theorem on the optimality of the Bayes classifier in this course; it is proved in the course of mathematical statistics and is not at all difficult. It says: if we knew the prior probabilities and the likelihood functions of the classes, that is, the class densities, we would write out the method of minimising the average risk, this functional, and the theorem states that the classifier for this must be built like this, as the $\arg \max$ of this expression. What we realised intuitively, on our fingers, is in fact the theorem on the optimality of the Bayes classifier.

The only thing we did not talk about then is ${\lambda }_y$, the losses from errors. We minimise the average risk, the expectation of the losses, where the size of the loss is ${\lambda }_y$, and it enters the optimal Bayes classifier as a factor. This is the only thing we could not have guessed intuitively. If you introduce a cost of error on the objects of the class $y$ and this cost, for example, has doubled, it is the same as if the probability of objects of this class appearing had doubled. That is, ${\lambda }_y$ and $P(y)$ form a product; ${\lambda }_y$ is in some sense a recalculated estimate of how often objects of this class will appear for you; they flow into each other. That is why in Bayesian classification theory the product ${\lambda }_{y}P(y)$ is always considered from here on.

And one more sardonic remark: of course it is not optimal at all; "the optimal Bayes classifier" is only a name. It is optimal only if the class densities are known exactly. As soon as we start estimating them from the sample, there can be no talk of any optimality: sample estimates of probability densities may well be overfitted. The Bayes classifier with estimates substituted into it is simply a way of constructing classifiers on the basis of a probabilistic model of the classes, nothing more.

And it seems that this is a bit of using a sledgehammer to crack a nut, because estimating densities is, generally speaking, a harder problem than classification. Why? Densities are estimated in a function space, a density is a continuous function, while in a classification problem there is a finite number of classes; classification is a piecewise-constant function, and recovering it seems simpler. Although both are recovering a function from data, so perhaps one cannot judge unambiguously that we are reducing the simpler classification problem to solving the harder problem of density estimation. People sometimes say so, but it is an interpretation, not a strict mathematical fact.

A survey of Bayesian classifiers

I will skip the piece about the naive Bayes classifier, although it is of course interesting, and go straight to the survey. We have to substitute density estimates into the construction of the optimal Bayes classifier, and today we have learned three types of them: parametric, nonparametric and mixtures of distributions. Let us do it.

The Parzen window method

Let us begin with the simplest case. If we substitute the nonparametric density estimate into the optimal classifier, we unexpectedly discover the Parzen window method, which we have already covered, a metric classifier; that is just how we wrote it out. The only thing we did not have then is the factor $P(y)/{\ell }_y$. It may cancel if the number of objects in a class is proportional to the prior probability of the class, that is, if there was no special distortion when the data were formed. But if, for example, the classes were unbalanced and you balanced them, evened them out, then this balance must be taken into account, in exactly this way. In the metric classifier there was no talk of this, so we did not know that such a factor had to be put in there too. It turns out that the Parzen window method, which we derived as a heuristic generalising the nearest neighbour method, is in fact statistically justified: it is the optimal Bayes classifier with the Parzen–Rosenblatt density estimate substituted into it. Such is the connection that has formed between two approaches that at first sight are about completely different things.

The quadratic discriminant

What if we apply the parametric approach, if the class densities are Gaussian? This leads to an interesting construction called discriminant analysis. It was studied a great deal from about the thirties to the sixties. Nowadays discriminant analysis is practically not used, because logistic regression and SVM certainly work better, but for general education it is good to know what it is.

An interesting point: if the densities are Gaussian, what will the separating surface look like? It is the locus of points where the expressions under the $\arg \max$ for two classes are equal. If we take logarithms, we get the first statement: the separating surface is quadratic. The quadratic form $x^{\mathsf{T}}\cdot x$ enters only in the construction with the inverse covariance matrix. If the covariance matrices of the Gaussian classes are the same, and that means the classes are of the same shape, their scatter ellipsoids are, as it were, moved apart by a parallel translation, identically oriented and with identical variances, then the quadratic terms cancel, and the separating surface degenerates into a linear one. This construction is called the linear discriminant, Fisher's linear discriminant: a linear separating surface, even if the classes are not spherical, as long as they are of the same shape, with the same covariance matrices.

On the slide there are six pictures with the geometric meaning of the quadratic discriminant. When we play with the factor ${\lambda }_{y}P(y)$, the separating surface simply shifts in parallel. The rule is easy to remember: it moves away from the most dangerous class, the one that either appears more often or has a larger $\lambda$, the loss from an error on objects of this class. If the classes are of different shapes, the surface is quadratic, and it can be so quadratic that a strange situation arises: a point that lies far away, closer to the objects of the other class, you will suddenly assign to the first class. This is an artefact of the method, an artefact of the generative approach: you estimate densities and then use them for classification, and far from the boundary between the classes anything can happen; the density estimate there may simply be inadequate, and for both classes at that. And situations where the classes are dissolved in each other are another matter: nothing can be done there; that is an unsolvable classification problem.

Fisher's linear discriminant

A little more about Fisher's linear discriminant. When you have described all the classes by densities with one and the same covariance matrix, you determine the estimate of this matrix, by and large, from the whole sample: the shape of the classes is the same, so it should be estimated from the full sample, only you shift each object $x_i$ by the expectation of its class $y_i$. Such is the small correction to how the covariance matrix then has to be estimated. It is as if you determined the centre for each class, brought the centres together at one point by a parallel translation, so that all the objects have shifted, and estimated the common covariance matrix from this sample.

The classifier comes out linear. To estimate it, one has to be able not only to compute the covariance matrix from the data but also to compute its inverse. Then by these formulas we compute the free terms ${\beta }_y$ and the coefficients of the linear classifier ${\alpha }_y$. For each class we get a linear function of its own, and the direction vector ${\alpha }_y$ of the $y$-th class says in which direction the objects of this class lie: the farther in that direction, the more numerous and denser the objects of this class will be. Here is one more geometric interpretation: for each class you have indicated the direction in which its objects are predominantly located. Clearly, the thing is not ideal. But it is one more example of a linear classifier. How many do we have already: SVM is one, logistic regression two, Fisher's linear discriminant three, and if you scrape around, a few more will be found in our course; many times we have recalled linear classifiers.

Gaussian mixtures and the radial basis function network

The last story I wanted to tell. Two stories I actually did not get to: about the naive Bayes classifier and about the Gaussian mixture with diagonal covariance matrices. Let us describe each class in the optimal Bayes classifier by a Gaussian mixture; I will show right away where we will end up.

We take all these estimates, which of course are obtained by the EM algorithm, substitute them, and unexpectedly discover an interesting construction called the radial basis function network, a three-layer network. The first layer computes the estimates of the mixture components for each class, the second layer mixes the components, and the third is the $\arg \max$, the optimal Bayes classifier. The density of each class is described by components, and each class may have its own number of components $k_y$. In architecture this is RBF, radial basis functions, and such a network is trained by the EM algorithm: for each class separately you recover a mixture of distributions by the EM algorithm. The number of components can be chosen there, which is the number of neurons in the first layer for each class, roughly as I described: increase it a little, decrease it a little.

Where this is applied

Is this applied now? Strangely enough, it is. Mixtures of distributions arise often and in many places. Should mixtures of distributions be used to recover classes in a classification problem? Sometimes they should: when we want to understand the internal structure of a class. There is a more general approach that really is applied: solving a classification problem through a clustering problem. First we understand the internal structure of each class, and this allows the classification problem to be solved more accurately. Gaussian mixtures for the optimal Bayes classifier are sometimes used precisely when something more than just classification is needed.

The second case is that very naive Bayes classifier that I did not have time to tell about. It is often used as a whipping boy, as a baseline: very simple, very easy to train, but it is a feeble classifier, based on the draconian assumption that the features are independent, so it rarely works anywhere. But sometimes it does work, for example in text analysis, where the features are word frequencies: there are very many words, and to a first approximation they can be considered independent of one another, so there the assumption can be stretched to fit, more or less.

In general, Bayesian classification theory is an achievement of somewhere around the sixties and seventies, and it still works in many places, wherever we need a generative approach. It may be needed for various reasons: it is not necessarily a strong classifier, but it does recognise the internal structure of the classes. Sometimes one needs to split a class into components and also estimate the probabilities of an object belonging to these components, that is, the reason why this object ended up in this class. Sometimes such things are needed for estimating risks, for obtaining structural information about what kind of object this was.

Conclusions

Seminar