But rather than start with the statement of Bayes’ Theorem, I want to use an old math teacher trick (which I realize many students hate) of trying to derive it from scratch, without stating what we’re trying to derive. Rather, we’ll start by modifying a problem that I described in an earlier post on probability distributions.

Lets pretend we have a robot arm with two joints: The first is fixed to the center of a table and spins horizontally. There’s a bar from this first joint to the second joint, which also spins horizontally and is attached to a second bar. The two bars are the same length, as shown on the left in the Figure below. The game is to randomly pick a pair of angles for the two joints, then try to guess the x and y coordinates of the hand at the end of the arm. As I described in the earlier post, we can think of the (density function of the) probability distribution of all the possible (x, y) coordinates, and it would look something like what’s shown on the right of the Figure. Here, darker colors indicate larger values of the function.

In the earlier post, we determined that the best place to predict that the hand would be is near the center, since that’s where the probability density is highest. But now we’re going to modify the game a bit: What if each time after spinning the wheel, we are told something about the x-coordinate; either its exact value or a small range. Then how would that change our prediction? Note that if we’re given the exact x-value and asked to predict y, then this is essentially the regression problem except that the probability distribution doesn’t look anything like the one we used in the post on regression.

But lets start with the case where we’re given a range, for example if we knew that the x-value was between 1/4 and 1/3. Then the probability of getting any point with an x-value outside this range would be zero, so the probability density that we would use to guess the y-value would be equal to zero for all points with x-values outside this range. So it would look something like the Figure on the right.

But dropping those values to zero isn’t enough to get the new distribution; the problem is that when we add this restriction on x, the probability of any point with an x-value in the correct range will increase. The question is: By how much will they increase?

In order to answer this question, we need to look more closely at what the probability density function really means. The first thing to note is that the value of the probability density function at a point is not the probability of choosing that point. In fact, the probability of picking any one point is zero, since there are infinitely many possible x and y values.

In order to understand the meaning of the probability density function, we need to use integrals, but (as usual) we can avoid much of the technical details by describing things in terms of the geometry that underlies those integrals. In particular, we’re going to think of our probability density function as describing the elevations of a mountain whose base is the square in which our robot arm rotates. But it won’t be a normal looking mountain – because of the way the density function looks, it’ll have a high peak in the middle, surrounded by a deep moat, then a high circular ridge (shorter than the central peak) around the outside.

I’ve attempted to draw this on the left, but you’re probably better off using your imagination to picture it. Once we’ve transformed our density function into this mountain, we can replace the word “integral” with “volume” and we’ll be able to calculate some probabilities.

Now, as noted above, if we pick one specific point, the probability that the hand will end up there is zero. However, if we pick a particular region A of the square, such as a rectangle defined by a range of x-values and a range of y-values, then there may be a non-zero probability that the hand will stop within A. (Though everything I will say below also holds true for more complex shapes A, as well as for shapes in higher-dimensional probability spaces.) In particular, the density function is defined specifically so that the probability will be equal to the volume of the part of the mountain above the shape A.

In other words, if we were to take a band saw to the mountain, following the outline of A, then the volume of the piece that we cut out would be equal to the probability of the robot’s hand stopping within that region. We’ll call this volume/probability P(A).

So, not only is the value of the probability density function at a point not the probability of getting that point (since it’s always zero), the value of the density function at a point doesn’t even need to be less than one. In particular, if there is a region A with a small area but a very high (though still less than 1) probability, the values of the density function would need to be very high in order to get the appropriate volume.

If we choose region A to be the entire square then P(A) = 1 because the arm is constrained to stay within that region. So the volume of the entire mountain is 1. If we choose a smaller region A, and an even smaller region B contained in A, then we’ll get a piece with smaller volume P(B) < P(A), and thus lower probability as we would expect.

But now lets return to the original question of how to modify our density function after we’ve narrowed down the set of possible outcomes to a smaller range A of x-values. This function will define a different mountain that is low and flat everywhere outside of A, but has the same elevations as the original within A, as on the left in the Figure below. We want to modify this function to give us a new probability density defining a volume function which we’ll write P(*|A) where * can be any region of the square. Since we know the robot hand landed in A, the overall probability, i.e. the volume P(A|A) of the new mountain, should be 1. However, since all we did was flatten the parts of the mountain outside the region, its volume is initially quite a bit less than 1.

In order to get the correct volume, we’ll need to scale the function up, i.e. multiply each value of the function by a constant k. The resulting function will define a mountain more like the one shown on the right above. For any region B contained in A, we’ll have P(B|A) = kP(B). Since we want P(A|A) = 1 = P(A)/P(A), the only possible value for k is 1/P(A) and we get P(B|A) = P(B)/P(A).

This is the case when the region B is contained in A. But what if it isn’t? For example, if we want to predict a range of y-values once we know the robot hand is in a certain range of x-values, then A will be a vertical strip of the square, and B will be a horizontal strip of the square, with the two intersecting in a smaller rectangle. (But as I noted above, we could just as easily let A and B be arbitrary blobs in the square or even blobs in a higher-dimensional space, but lets not get too complicated…)

So, if we want to calculate P(B|A) in this case, we need to consider two parts of B separately: The density function above the part of B that is outside of A will all get flattened to zero, so P(B|A) is completely determined by the part of B inside of A, i.e. the intersection A ∩ B. In other words, we have P(B|A) = P(A ∩ B)/P(A). Note that there’s no difference between A and B in this formulation, so we also have P(A|B) = P(A ∩ B)/P(B). We can solve both equations for P(A ∩ B) to get P(B|A)P(B) = P(A ∩ B) = P(A|B)P(A). Finally, if we divide both sides by P(B), we get Bayes’ Theorem:

P(B|A) = P(A|B)P(A)/P(B)

Of course, for the problem we started out with, the original equation P(B|A) = P(A ∩ B)/P(B) may sometimes be more useful. But in the standard setting of Bayes’ Theorem, P(A ∩ B) is the probability that both events happen (or both statements are true) so it might be harder to calculate.

For extra credit, take a minute to think about how you might calculate the probabilities of different y-values if we knew the exact value of x rather than a range. I’ll give you two hints: First, note that the probability density function over the vertical line defined by a single x-value defines a single-variable function like you might find in Calculus I and II, and there is some area (rather than a volume) below this function. Second, note that you can take A to be a small rectangular strip around the line defined by the x-value, calculate its volume, then make the strip smaller and smaller and take a limit.

But this post is already long enough, and I expect that most of my readers either don’t want to read about limits, or would rather work it out themselves. (For me, it’s both.) So I’ll leave it there.

The goal of symbolic regression is to optimize a non-linear regression model. The outcome of the process will be an equation that defines a function from our input features to an output value. For example, if we have two input features –  and – then the model may look something like . To measure the quality of a given model, we would use it to predict values for a training set, then add up the squares of the differences between the predicted values and the actual values. We want a model for which this sum is as small as possible.

Notice that the plus signs in the equation naturally decompose it into pieces. In my earlier post on non-linear regression, we approached this problem by deciding what kinds of pieces we wanted to consider – for example, all polynomials with powers up to 3 – and then using standard optimization (gradient descent, though I didn’t call it that at the time) to select the best coefficients. This poses the risk of getting stuck in a locally optimal solution, but there’s also another problem: When we add these different types of elements into the equations, we’re effectively creating a higher-dimensional problem. If we want to allow a lot of different type of elements in the final equation, things can become computationally unruly. For example, if we start with 10 features instead of 2, then there are 59,048 possible polynomials with powers up to 2. If our data has 100 features then the number is roughly a 5 with 47 zeros behind it.

But the subtle point here is that even though we want as many options as possible for the types of elements allowed in our solution, we’ll usually want the solution itself to contain only a small number of elements from this large set of possibilities. This, as it turns out, is also something that a genetic algorithm excels at.

As the name suggests, the idea behind a genetic algorithm is to mimic the natural selection process that produces genetic codes, i.e. DNA. It’s common to think of DNA as a collection of genes whose presence or absence determines the traits of a plant or animal. This is a gross oversimplification of the biology, but that’s never stopped us before.

To solve an optimization problem with a genetic algorithm, you start by defining the genes of the thing that you want to optimize. For symbolic regression, the “genes” are simply the pieces of the equation between the plus signs. So the equation from above has four genes – , , and .

In the biological world, we get new genetic codes when existing codes combine with each other, mutating slightly in the process. We will once again ignore the biological complexity and define our own simple ways to combine and mutate the things that we’re interested in. There are a few different ways we could combine two equations, such as randomly picking a few elements from each of the original equations and combining them. So if we start with equations and , we might combine them to form .

Next we need a way to “mutate” our equations, and again there are some fairly simple options. For example we could randomly delete and/or add a small number of elements to the equation, or we might adjust one or more of the coefficients slightly. The combined equation might become . Note that even if each equation that we start with only has a small number of elements, we can randomly choose elements to add from an arbitrarily large set of options without running into computational problems.

Once we know how to combine and mutate our regression models, we can start the Darwinian process. We’ll begin by randomly creating a small initial “gene pool” of regression models. Then we’ll take a pair of models, combine them, mutate the result, add it back into the original gene pool and repeat.

The question, of course, is how to decide which pairs of models to combine each time. Survival of the fittest implies that only the “best” models should be allowed to pass on their genes, and as described above, we have a way to determine how “fit” each model in our gene pool is. However, if we always combined the best two models, we would likely fall into the very trap that we were trying to avoid; even though there would be some amount of randomness from the combination and mutation, for the most part we would get lots of models similar to the top two of the models we started with. These might all be centered around a locally optimal solution that is not globally optimal.

So to fix this issue we’ll add in a third element of randomness: Rather than always selecting the best two models to combine, we will randomly pick them, but with a bias that makes the better models more likely to get selected. That way the models in the gene pool will tend towards better solutions, but they will also spread out randomly and explore possible solutions that are quite a bit different from the best solutions.

An additional advantage of this approach is that it can be done in parallel; while gradient descent always looks at a single variant of the one best-so-far solution, a genetic algorithm can explore an unlimited number of different possible solutions all at the same time.

Genetic algorithms are commonly used for feature selection, where each “gene” is a feature in a high-dimensional data set, and the goal is to select a small set of features determining the best lower-dimensional representation of the data set. For this kind of application, the tough part is deciding what measure to use for the “fitness” of each set of features. I used symbolic regression as my example above because the fitness function is straightforward. Though, in fact, you can think of symbolic regression as a form of feature selection in a high-dimensional kernel space for a low-dimensional data set.

Of course, there are some optimization problems that genetic algorithms can’t help you with, particularly if you can’t break your problem down into “genes” such that a solution is uniquely determined by which genes are on or off. For example, if you’re training a complex mixture model on low-dimensional data set, a genetic algorithm probably won’t help you much. You can deal with this to some extent by running gradient descent many different times from different random starting points. (This also has the advantage that you can run these all in parallel.) There are also more general randomization techniques such as simulated annealing that don’t require a notion of “genes”, though they do have other constraints. But those will have to wait for a future post.

Lets say that we have a horizontal robot arm that is attached to a square table by a joint that can spin 360 degrees. The arm consists of two pieces, connected by a second joint that can also spin horizontally 360 degrees, as on the left in the Figure below.

Here’s the game: This arm is behind a screen so that we can’t see it. The arm is programmed to randomly select two angles, then move its joints to those angles. Once it’s done, it beeps and we have to guess the coordinates of the hand at the end of the arm, with respect to the table. This is, of course, going to be difficult since we have essentially no information to go on. However, we’re not completely in the dark because we know from the mechanics of the arm that there are certain places it is more likely to be than others.

If you ever played with spirograph (a toy from when I was a child that I was pleasantly surprised to discover you can still buy today) you know that when you move in two circles like this, you spend much more time near the outside and near the center than in the area in between. So the probability distribution defined by the possible locations of the hand looks something like the picture on the right of the image above: It is darker along the edge of the circle of possible locations of the hand, as well as in the very center of the circle, both areas where the hand is more likely to be. So if we wanted to guess the coordinates, our best bet would be right in the middle of the table.

As I hinted at above, this robot arm is an example of what’s called a configuration space: We have a number of variables (the angles of the arms) that define different possible configurations. The configuration space is a way of thinking about all possible configurations of a system (such as the robot arm) in a unified way.

This is similar to the idea behind the data spaces that we’ve seen in past posts, made up of all possible data points, but with one important difference: The data spaces we’ve seen were all vector spaces, with their points uniquely defined by their coordinates. In this configuration space, however, each configuration is defined by two angles. The set of all possible angles make up a circle rather than a line, so the “shape” of the configuration space is actually a two-dimensional torus (a doughnut) like the one shown to the right. You can think of this as the result of dragging the smaller circle defined by the elbow angle around the larger circle defined by the center angle.

Meanwhile, each configuration in the configuration space defines a vector (given by two coordinates) in the vector space defined by the table. Note that except for the points right at the inner edge and the outer edge, each vector is actually defined by two different configurations/positions of the arm: one in which the elbow joint bends to the left and one in which it bends to the right.

In topology, this is called a map or function: a way of identifying each point in one space (the doughnut) with a single point in some other space (the vector space). It is similar to the traditional notion of a map as a piece of paper in which each point in some geographical region is identified with some point on the paper. Just as paper maps sometimes distort the relationships between the points that they’re representing, such as when one makes a flat map of the earth, the map defined by the robot arm distorts different parts of the doughnut-shaped configuration space in order to fit it into a flat square. This distortion is recorded by the probability distribution, which becomes denser in the places that are squeezed together (near the outside and inside of the region) and less dense in the places that get stretched.

But what does this have to do with the sorts of data analysis that I’ve described in the past on this blog? If we’re trying to predict what restaurant someone will want to go to for lunch or whether or not a text message is spam, there isn’t a mechanical configuration space behind a curtain determining the answer.

However, from a more philosophical point of view, it can often be useful to think of situations in this way. For example, let’s imagine hypothetically that there was a configuration space of all possible people. This certainly wouldn’t be a vector space since people are much more complicated than that. And it would be infinitely more complex than the doughnut defined by the robot arm. But the actual structure of this configuration space doesn’t matter. What does matter is that when we start recording data about people, such as their food preferences, we are essentially building a map from this (hypothetical) configuration space to the data space.

Because we don’t know the structure of this “configuration space of people” the way we understand the configuration space of the robot arm, we can’t explicitly determine what the probability distribution defined by the way this map stretches and squishes points looks like. That’s why we have to build an approximation of it based on that data points that we think of as being sampled from it.

But the key to this perspective is that when we have a map from a configuration space to a data space, a lot points in the configuration space get squished together in the data space. This happens with the robot arm because each vector of coordinates comes from two different arm positions. So, for example, if we knew the coordinates and wanted to predict the angle of of the elbow, there would be two possibilities, each of them equally likely. In other words, given our knowledge of the coordinates, there would be a 50% chance of either of the two possible angles being correct.

So this gives us one way to interpret probabilistic statements about the present. But what about the future? Well, if we’re willing to believe that a given configuration space is deterministic, i.e. that the configuration in ten minutes is completely determined by the current configuration, then a probabilistic statement about the current configuration can easily be translated into a probabilistic statement about the future. For example, if we knew that each joint in the robot arm rotated at a constant (known) speed then knowing its current configuration would tell us what the configuration would be in ten minutes. But if all we know is the coordinates, then there are two possible configurations that it could be in in ten minutes. So that would give us a probabilistic statement about the future.

A better example of this is weather forecasting. (And it was while reading the weather chapter of The Signal and the Noise that I really started thinking about this.) The weather is determined by the movement of massive numbers of molecules in the Earth’s atmosphere, by photons arriving from the sun, and a handful of other physical phenomena. From the perspective of Newtonian physics, this makes the weather one giant configuration space – Just like the robot arm, but instead of the moving parts being two joints, the moving parts are all of the particles in the atmosphere bouncing off of each other.

According to Newtonian physics, this model is deterministic: If we could know the positions and velocities of all these particles at a given time, we could (in theory) determine exactly what their positions and velocities will be in ten minutes, or better yet tomorrow morning. Quantum physics, on the other hand, suggests a bit on non-determinism, but lets not worry about that right now.

When weather forecasters create their models, they can’t keep track of every single air molecule. Instead they take measurements from a relatively small number of weather stations throughout the world. These measurements define a map from this massive and unwieldy configuration space to a data space defined by the temperature readings, etc.

So, when NOAA states that there is a 20% chance of rain next Tuesday, they are in effect, saying the following: Of all the possible configurations that the atmosphere could be in, given these readings, 20% of them would lead to rain on Tuesday and the other 80% wouldn’t. In other words, the there is a 20% chance of rain, given what is known about the weather configuration space. In theory at least, that’s what we should expect other types of forecasts to mean as well, though as The Signal and the Noise beautifully explains, in practice some types of forecasts are better than others.

The basic idea behind the K-means algorithm is essentially a very restricted mixture model: We select a number K and create a mixture model with K simple models, each of which is a Gaussian blob. However, unlike in general mixture models where we allow each Gaussian blob to be a different size and to be stretched in different directions, as on the left in the Figure below, K-means requires that each blob be a fixed size and completely symmetrical (so no stretching), as on the right. You can imagine the sorts of complaints that this would face: Most data sets aren’t made up of symmetric Gaussian blobs of the same size. In fact, unless you miraculously scale your data perfectly, you’re pretty much guaranteed that it won’t be.

However, if the data is made up of well separated clumps, even if those clumps aren’t perfectly round and all the same size, K-means will often (though not always) be able to find them. (These “clumps” are often called clusters, but I’m going to leave a general discussion of clusters for a later post.)

So when you use K-means as a substitute for a mixture, you’re taking a risk that the final answer will be pretty lousy. The trade-off for taking this risk is that the algorithm is much much faster than training a general Gaussian mixture model. K-means still follows the two-step Expectation Maximization outline of a mixture model, but each step has been optimized to the point that you can barely tell where it came from.

First, note that since each of the Gaussian blobs in our mixture model is the same exact size and shape, the only information we need to record about each blob is its center. So our “model” will consist of K points spread throughout the data space. Recall that the first step in the Expectation Maximization algorithm is to assign each data point to one of the simple models (i.e. one of the Gaussian blobs) based on which one assigns it the highest probability. The probability that one of the Gaussian blobs assigns to a given point goes down as we get farther away, and all the blobs are identical. Thus each data point will simply be assigned to whichever center point is closest to it. That’s step one.

To understand what’s happening here, consider the data points shown in black in the figure below. We have a smaller number of blob centers, shown in blue. When we assign each data point to its closest blob center, this should remind you of the nearest neighbor algorithm. In the post on nearest neighbors, we noted that the set of all points that were closer to a given data point than to any other data point form what’s called a Voronoi cell around each data point. For K-means, we’ll get something very similar, except that the Voronoi cells are around the blob centers rather than data points. The Voronoi cells are outlined in blue on the left of the Figure. Each data point will be in the Voronoi cell around one of the blob centers, and we will assign it to this blob.

For step two, we need to adjust each of the Gaussian blobs to maximize the probability of all the points assigned to it. In general mixture models, this involves repeatedly adjusting the parameters using the notion of a gradient. However, if we’re not allowed to stretch the Gaussian blobs then we can immediately figure out exactly where to put them: We’ll want the point where we get the smallest possible number when we add up the distances squared from the center point to each of the data points assigned to it, i.e. all the data points in a given Voronoi cell.

In other words, if we take all the data points in one of the Voronoi cells and think of these data points as heavy balls connected by bars that hold them in place, we will want to move the center point of the Gaussian blob to the center of gravity of this structure. This center of gravity is called the centroid. As it turns out, the center of gravity is straightforward to calculate – you just take the average of the data values in each dimension. So for step two we move each Gaussian blob so that it is centered at the centroid of the points that were assigned to it in the first step. These centroids are shown in red in the center Figure.

Then we repeat – We go back through the data points and reassign them based on the new center points of the Gaussian blobs. (In other words, we form new Voronoi cells around the red points.) Then we calculate the centroids of each of the new subsets of data points (the points in each new Voronoi cell), then reassign the data points and so on.

After we repeat these two steps enough times, we’ll get to a point where step one assigns the data points to centers in exactly the same ways as the previous time that we ran it. So there won’t be a need to calculate the centroids again (since we’ll just get the ones we already have) and the K-means algorithm will terminate. In the Figure above, these would be something like the blue points shown on the right. It’s also common to run through the two steps a fixed number of times, rather than running it all the way to the point where the centroids stop moving, hoping that this will give you a reasonably good approximation of the final answer.

If the data happens to be separated into a certain number of well defined, tightly packed clumps/clusters and you happen to select K to be equal to this number then there’s a reasonable chance that the blob centers that the algorithm ends with will be the centers of these clumps/clusters. As always, we usually won’t be able to directly “see” whether or not this is the case, since most data has more than three dimensions. However, there are a number of measures that have been devised to measure how good K-means has done.

The most common one is to compare the distances between the final center points to the average distance from each data point to the center that it is assigned to. If the first of these numbers is much higher than the second then that means that K-means has found well separated clusters. Since the quality of the final answer depends heavily on picking the right value of K, this is often used to decide which K to use: Simply run K-means with different values of K and choose the one in which the difference between the two numbers is as large as possible.

Of course, this doesn’t solve the problem of what to do when our data set is not composed of a small number of well separated clumps. One of the main criticisms of K-means is that many (perhaps most) interesting data sets tend to have more complicated structures, so it would be naive to try using K-means on them. But like I said at the beginning of this post, I’ve now come to like K-means, nonetheless, by looking at it from a different perspective. In particular, I want to suggest using K-means for an off-label usage very similar to a technique from signal processing called vector quantization. This use of K-means has been suggested by others in the past, but doesn’t seem to be widespread. (It seems to be better known within the signal processing community, since that’s where the idea of vector quantization is well understood, but not as much within machine learning and data mining.)

Notice that the data set shown in the Figure above is not made up of a small number of round Gaussian blobs. In fact, it isn’t even made up of stretched Gaussian blobs, since one of clusters/clumps of points forms a curved shape. But you can see that the final blob centers found by K-means still reflect the structure of the data set, just in a more subtle way: We end up with three blob centers within the round cluster in the lower half of the data set and the rest of the blob centers forming an arc that follows the curved cluster.

To understand why this happens, note that each step in the K-means algorithm reduces the average distance squared from the data points to the blob centers. In other words, we can think of K-means as trying to place the blob centers so that each data point is as close as possible to some blob center. This has the effect of spreading the blob centers evenly around the data set, so the “shape” of the set of blob centers will be very similar to the overall shape of the data set.

So the blob centers found by K-means are very similar to a random sample of the original data set. However, there are two important differences: First, unlike a random sample, the blob centers are not in the original data set. One slight modification would be to take a sample of the original data set by choosing the data point closest to each blob center. Second, the blob centers are not random (though they do depend on the initial choice of centers, which are often chosen randomly.)

The blob centers have two important advantages over the original data: First, there are fewer of them, allowing one to run more sophisticated algorithms in a reasonable amount of time. Since K-means on its own is very fast, this would be a good trade-off. Second, the set of center points are less random than the original data set. In data analysis it’s common to think of randomness as being a good thing: We want random samples of our data to make sure it’s evenly distributed, and of course the field of Statistics is devoted to understanding all the facets of randomness. However, from geometric perspective the problem with random samples is that by their nature, they contain large gaps that can be mistaken for geometric features. This is a fundamental principle of randomness, and in fact you should be suspicious of any “random sample” that does not have both many large gaps between points and many points that are very close together.

One idea that has come out of numerical analysis (a field within mathematics) is that a homogeneous sample is often better than a random sample for large scale simulations. In other words, it’s often good to have data points that are “evenly” spread out within the probability distribution, but with more concentration in the darker parts of the density function. If we choose the number K of blob centers to be relatively large (though still much smaller than the number of data points) then the final set of points that K-means gives us should have this property: It will be relatively evenly spread out, but still have the rough shape of the original data set. So the set of blob centers can be used for the types of more computationally expensive forms of exploratory data analysis (which we’ll discuss in upcoming posts) whose goal is to understand the overall structure of the data well enough to decide what types of models and what types of kernels to use.

Another possible use for the blob centers is to use them as the anchor points for a Gaussian kernel. Recall from last week’s post that a Gaussian kernel is defined by choosing a small number of points in the data space. We want these points to be fairly evenly spread out, but close to the original data set. If we use all the original data points for make the Gaussian kernel then we are more or less guaranteed that our model will suffer from over-fitting. Thus the blob centers will often be a much better option for defining a Gaussian kernel.

So, that’s my take in K-means: If you think of it as a clustering algorithm, or as a substitute for a Gaussian mixture model, then you may be asking for trouble. However, I think that K-means has a great deal of unexplored potential as a quick way to reduce the data set to a relatively small number of representatives for us in exploratory data analysis or forming a Gaussian kernel.

The idea behind a mixture model is that the underlying distribution for a data set may have different types of local structure in different parts of the data space, as in the Figure below. For example, if different data points are the result of different phenomena (for example, different species in a data set of animals, or different demographics of customers in a marketing data set) then we would expect each phenomenon to be described by a different model. We can get a model for the entire data set by combining the models for all the species/demographics/etc.

We’ve now seen a number of different models that we might want to combine. There’s the linear regression model, where we choose a line/plane/hyperplane and define the density/probability of a point to be the Gaussian function of the distance from the point to the line/plane hyperplane. In other words, the probability is highest right on the line/plane/hyperplane and gets closer to zero as you get farther away. Such a distribution is shown cutting diagonally across the Figure above. We later saw the logistic distribution, where we again choose a line/plane/hyperplane, but this time the density/probability gets closer to one as you get farther away on one side of the line/plane/hyperplane, but gets closer to zero as you get farther away on the other side.

A third type of model/distribution that we saw, but I didn’t describe in quite as much detail, is the Gaussian blob.  The basic Gaussian blob is what you get by choosing a point in your data space and defining the densisty/probability of the distribution to be the Gaussian function of the distance from the point. So, the density is highest right at the point, then gets closer to zero as you get farther away. A Gaussian blob is similar to the regression distribution, except that it’s based on a point rather than a line/plane/hyperplane. You can then create more sophisticated Gaussian blobs by stretching this blob in different directions. The Figure above shows two Gaussian blobs, on on either side of the regression distribution. This was the distribution that was implicit in the PCA algorithm from a few posts back.

We can also take any one of these three models and apply a kernel to it. In fact, we could apply a different kernel to each individual simple distribution in a mixture model before mixing them, and/or we could apply a single kernel to the entire mixture at the end of the process. However, lets not make things too complicated for now.

The next step is to combine these different distributions,  and there are two basic ways we can do it. Recall that each simple distribution defines a probability/density value for each point in the data space. We need to use these different density values at each given point to decide what its probability/density in the mixture model should be. The first possible way is to add together the values the densities of the simple distributions. This is roughly what we would get if we printed each distribution on a sheet of clear plastic (like an old-fashioned overhead projector slide) and then stacked them all on top of each other: The overlaps between the densities add together to give us the overall density.

Another way to think about this is to think of the simple distributions as piles of sand. From this perspective, a regression distribution looks like a long, narrow ridge, something like what a gopher makes when it digs a tunnel near the surface of your lawn. A logistic function will look like the edge of a plateau – very flat on either side, but with a steep wall going up from the lower side to the higher side. A Gaussian blob will look like a single mountain peak, though its base may be an ellipse rather than a perfect circle. When we add the distributions, we pile them on top of each other like the middle picture on the Figure below, where we combine the two Gaussian distributions shown at the top of the Figure.

The second possibility is to take the maximum density of all the simple models. To picture this, we’ll again think of each simple distribution as a mountain whose height is determined by its density. To combine by the maximum method, we’ll place them all in a mountain range, as if they were holograms that can pass through each other rather than piling on each other, as in the bottom picture in the Figure above. With this method, each simple distribution will peak above the others in certain places, but will be below other distributions in other spots. The density of the mixture distribution will be the overall height of this range at any given point.

These two ways of combining distributions seem very different in how they treat the overlaps between distributions. The addition method makes more sense from a theoretical perspective: If each distribution is a separate phenomenon that is producing data points, then each simple distribution has some chance of putting a data point in any given spot, so the probability that we’ll get a new data point in that spot should be the sum of all the simple distributions. However, (as we’ll see below) the maximum method is much easier to deal with in practice because when we look at a single point, we can essentially ignore the effects of all but one of the distributions. Luckily, if the distributions that we’re combining are reasonably far from each other then their density values along the overlaps will be low enough that there often is not a huge difference between the two methods. So, it is common stick with the maximum method for ease of computation.

So, now that we know how to combine our simple distributions, we have two problems: The first is to decide how many simple distributions to combine and what types to pick. If you happen to already know something about the data set, such as a theoretical framework or some other form of domain knowledge, then this might be used to select the types of simple models. Otherwise, selecting models can be a difficult and subtle problem, similar to the problem of choosing the appropriate kernel for a data set. This problem generally falls under the category of unsupervised learning (the machine learning term) or descriptive analytics (the data mining term). In a few weeks, I’m going to switch gears and start covering unsupervised learning/descriptive analytics, but for now lets assume that we have somehow managed to choose the types of simple distributions we want to combine in our model.

The second problem is to choose the parameters for each simple distribution. Recall that the parameters for a regression or logistic distribution determine the line/plane/hyperplane that the distribution is based around. For a Gaussian blob, the parameters determine what directions the blob will be stretched in, and by how much. The set of parameters for the mixture model is the union of all these parameters for the different simple models. Moreover, if we run any of the distributions through a kernel, then that will give us even more parameters for the mixture model.

Recall that in the post on regression, we saw how we can fit a single model to a data set by trying to maximize the values of all the data points with respect to a probability density function. Specifically, we chose the parameters that give us the largest number when we multiply together the density values of all the points in the data set, giving us the distribution that is “most likely” to have produced the data set. This is what we called maximum likelihood estimation.

In practice, we often can’t calculate this value for every single possible distribution, so we’ll use what’s called gradient descent: We start with a randomly chosen distribution, then figure out what direction to nudge the parameters in to improve the probability score a little. (The gradient is what tells us how to do this.) Then we make these changes to the parameters and again figure out how to further nudge the new parameters of this new distribution, and so on. The process stops when we can no longer make small adjustments that improve the score. We saw this approach in the post on logistic regression, and it is very similar to the basic method for training neural networks.

We’ll select the parameters for a mixture model by a closely related algorithm called Expectation Maximization. We start by choosing initial parameters for each or our simple distributions. We can choose random parameters, though if we have another way of getting an initial estimate (such as domain knowledge) this will often give is a better final answer. Then we want to start adjusting these parameters. However, when we determine how to adjust the parameters, we don’t want to use all the data points for each simple distribution. From a theoretical perspective, most of the data points will have come from other simple distributions in the mixture, so each simple distribution should only care about the data points that it is likely to have produced. From a practical perspective, if we use all the data points for each distribution, we’ll probably end up with all the simple distributions having the same parameters, rather than having them nicely spread throughout the data set.

So, how do we determine which distribution was mostly likely to have produced each data point? This choice is an example of what’s called a latent variable in statistics: The value is not part of the initial data set and is not part of the final answer, but we need to choose a value in order to determine the final answer. It’s sort of like the middle layers of neurons in a neural network, playing a major behind-the-scenes role. In fact, these middle neurons were essentially how neural networks solve the equivalent problem, by having the middle neurons decide how to adjust the parameters of each earlier neuron when looking at a new data point.

But for mixture models, we’ll use a much simpler method than what we saw with neural networks: For each point in our data set, we will calculate its density/probability in each of the simple models with our initial parameters, then assign it to the simple distribution in which it scores the highest. The Figure below shows how the we divide points in the data space between the three simple distributions in the mixture model from the beginning of the post. Points in the red and blue regions get assigned to the Gaussian blobs, while points in the green region go to the linear regression distribution.  Essentially, this method treats the data points as if they were produced by combining the simple distributions using the maximum method, since it only allows each data point to be attributed to one distribution. This is the first step in the Expectation Maximization algorithm.

For the second step, after a collection of data points has been assigned to each distribution, we have each simple model pretend that this collection of data points was the entire data set, and adjust its parameters accordingly. In other words, we figure out how to make each simple distribution better fit the points that are assigned to it  (i.e. with gradient descent) and make these adjustments.

Once we’ve done this for each of the simple distributions, we note that the way we divided the data set between the different simple distributions may no longer be up to date: Some of the distributions may have shifted so that a data point now gets a better probability/density value in a different simple distribution. So, we return to step one and re-divide the data points. Then we re-adjust the distributions accordingly, and so on. At some point, the way we divide the points between the simple distributions will become consistent and we will run out of ways to improve the parameters. At that point, we’ll stop and use whatever collection of parameters we have for the different simple distributions to define the overall distribution.

As with any of these algorithms, we should consider its benefits and limitations. First, note that unlike regression and classification algorithms, mixture models do not explicitly predict anything about new data points. Instead, they give a nice description of data sets with more complicated structure. In particular, because they have at their core a small number of relatively simple models/distributions, they can produce a description of a data set that could potentially be interpreted by a human. In many situations, this will be more valuable than a black box classification or regression algorithm.

There are two major technical limitations: First, the initial values that we choose for the parameters can make a big difference in the final result. Since each step in the expectation maximization procedure only makes incremental changes, it can get stuck in configurations in which there is no small change that will improve things, but there are are larger changes that would produce a much more accurate model. Think of a ball rolling down a hill that gets caught in a small divot half way down. In order to get it to roll the rest of the way down, you have to push it up slightly first. However, the Expectation Maximization algorithm has no procedure for nudging a set or parameters out of analogous divots.

The second limitation is that the quality of the final model/distribution is very dependent on the choice of the number and types of the simple distributions in the mixture. As I mentioned above, if you don’t have a good external reason to choose a particular collection of simple models, this can be very difficult to determine. This is one of the motivations for unsupervised learning/descriptive analytics, which I’ll begin discussing in two weeks, after a brief digression next week into Gaussian kernels.

One of the main goals of PCA is to tell the difference between the three data sets shown to the left. If you run linear regression on each of them, you’ll get the same line, shown in blue. But there’s clearly a difference in how well this line approximates each data set – It does a really good job of summarizing the data set on the right, but a pretty miserable job for the data set on the left. For the middle data set, it’s somewhere in between. We can see the difference for a two-dimensional data set, but as usual we won’t be able to visually inspect the data like this if there are more dimensions. That’s where PCA comes in.

PCA can tell the difference between these three sets because rather than fitting a line or a plane or a hyperplane to the data, it (roughly speaking) fits an elliptical Gaussian distribution, such as those shown in the second picture. This distribution is defined by its axes or principal components – the red arrows shown in each picture. The angles of these vectors tell you the angle that the cloud is tilted at and the lengths of the vectors tell you how much it is stretched in each direction. The arrows on the left are almost the same length, but on the right, one of the arrows is much longer than the other. This large difference in length is what tells us that the data set on the right is much closer to a line along the longer vector.

This is probably a fitting place to remind the reader of one common misconception: PCA and any other computational/statistical method only measures correlations in the data, not causation. So no matter how closely the data fits a given model, I would be wrong to conclude that one factor or set of factors causes another. If two things are correlated, it might that one causes the other or it might mean that both are caused by an external factor not represented in the data. But a correlation could also be a coincidence or the result of a bias in the way the data was collected. One can sometimes establish causation by running controlled experiments, but no amount data analysis alone can turn a correlation into a causation.

OK, now that that’s out of the way, lets consider higher dimensional PCA. Because the examples above are two-dimensional data sets, there are two axes and two vectors. For higher dimensional data sets, there will be more axes. In particular, for an n-dimensional data set there will be n principal components. Some of these will be very short and some will be long, and the interpretation of long and short will be the same as in the two-dimensional case. So even though we won’t be able to see a high dimensional data set the same way we see a two-dimensional data set, we can get a rough idea of its shape from the principal component vectors.

For example, if we have a fifteen-dimensional data set with three really long PCA vectors and twelve really short vectors, then that means the data is mostly concentrated along a  three-dimensional hyperplane in a fifteen-dimensional space. Moreover, those three long principal components define this three-dimensional hyperplane along which the data lies. As in the picture, each principal component usually won’t be along one of our original features (though occasionally they might.) Instead, they will define some interesting correlations between combinations of features. (But again, we should never conclude that the value of one feature determines the value of another just because they’re strongly correlated.)

Unlike the form of regression that we discussed previously, PCA does not work by formally fitting a model to the data. Instead, it works by calculating the principal components directly from the structure of the data. The PCA algorithm is usually described in terms of variances and covariances between the different features/dimensions that define the space that the data lives in. But I want to explain it from a different perspective that leads to the same algorithm, but is a more geometrically intuitive. It also involves ideas that will come up later on in this blog.

We start by centering the data in our n-dimensional space. In other words, we find the center of mass in each direction, then translate the entire data set so that the center of mass in each direction is at exactly zero.

Next, we want to think of each data point as pulling on the data set, attempting to stretch it away from the origin, as indicated on the right. (The green arrows are the axes.) Note that the data points aren’t pulling away from each other – they’re pulling away from the center of mass. The farther a point is from the center of mass, the harder it pulls. With each point pulling together, we want to consider the overall effect that they have on the set. In particular, if there are a lot of points far away from the origin in one direction, as in the data set on the right of the above figures, the data set will mostly get pulled in that direction. Because it’s a stretching action rather than a translation, points that are pulling in opposite directions are actually working with each other rather than against each other.

The forces of all this pulling can be summarized in a certain matrix (which happens to be equal to the covariance matrix from statistics.) A direction in which the points are pulling the hardest is called that first eigenvector of the matrix and this will be the longest principal component. Next we look at all the vectors that are orthogonal/perpendicular to the first eigenvector and choose a direction with the most force among these, which we call the second eigenvector. The we look at all the vectors that are orthogonal/perpendicular to the first two eigenvectors and so on until we run out of orthogonal vectors.

The length of each eigenvector will indicate how strongly the data is pulling things in that direction. In practice, these vectors are calculated by a process called diagonalizing the initial matrix. So the PCA algorithm consists of calculating this matrix, diagonalizing it, then translating the results into principal component vectors.

The PCA gives you a lot more information than linear regression, but note that it does not have the flexibility of the types of regression we discussed with curved hyperplanes. So if your data follows a curved shape rather than a hyperplane, PCA may give you lousy results. There are ways to fix this, such as something called a kernel, but that will have to wait for a later post.

From the perspective of regression, the most important thing about high dimensional data is that the more dimensions there are, the easier it is to find a hyperplane that is close to all the points and the less meaningful such a hyperplane is. We know that for any two points in a plane, there is a line that goes through both of them. Similarly for higher dimensions, if the number of data points is equal to the number of dimensions then there will always be a hyperplane that contains all the data points. Moreover, if there are fewer points than dimensions, there will be infinitely many such  hyperplanes. Keep in mind that when all the points lie in the hyperplane, any error in one of those points will be amplified as hyperplane gets farther from it. (Think of it like a lever.) This means that any answer given by regression on such a data set should be viewed with suspicion.

As it turns out, data sets where the number of data points is less than the number of dimensions are quite common. For example, one of the current (as I write this) Kaggle competitions has 283 data points in the training set and around 28,000 dimensions (7,000 features extracted from each of four handwriting samples from each person.)

The second problem with high dimensional data is that because of the way distance is calculated, there are fewer points near any given data point, and a lot of points of roughly the same distance. For example, the ratio between the volume of a (high dimensional) ball of radius one and a high dimensional ball of radius two very quickly gets larger and larger as the dimension increases. In other words, the volume in the ball gets more and more concentrated near its outer boundary. You can check this by working out the volume formula for higher dimensional balls, but I won’t make you – Instead, lets consider a more concrete example.

We’ll start with a fairly simple picture: Take a square and in each corner draw a circle so that the four circles are as large as possible, just barely touching each other. In the middle of these four circles, draw a fifth circle (shown in blue in the Figure below) so that it just touches the four corner circles. It’s not too hard to calculate the radius of this middle ball if you remember the two-dimensional distance formula, but lets not worry about that right now. (Feel free to work it out on your own.)

We can do something very similar in three dimensions: Take a cube and put a sphere in each corner so that they just barely touch each other, as in the middle picture below. There are eight corners, so there are eight spheres. If we look at the cube from the side, like on the right of the picture, it looks a lot like the original square, except that each circle has become two spheres, one behind the other. We can also put a sphere between the eight corner spheres and make it as large as possible so that it just barely touches them. You can’t see the blue sphere in the middle picture because it’s behind one of the corner spheres, but you can see it on the right.

Notice that in the picture on the right, the blue sphere appears to overlap the spheres in the corners. Even though it just barely touches them, the points where they touch are just behind the corner spheres. So when we project the three dimensional configuration into two dimensions, they appear to overlap. If the square and the cube both have sides of length one then the eight corner spheres will have the same radius as the four corner circles, namely . However, the blue sphere will have slightly higher radius than the blue circle. (Again, you can calculate the exact radius using the three-dimensional distance equation.) So the corner spheres in three dimensions are farther from the center than the corner circles in two dimensions.

What would happen if we did this in dimension four? In other words, we will take a four-dimensional cube and put sixteen (higher-dimensional) spheres in its corners so that they just barely touch, then place a blue (higher dimensional) sphere between them that is just big enough to touch the sixteen corner spheres. I don’t expect you to be able to picture this, but hopefully you can believe that if we lived in four-dimensional space this would be possible. In fact, we can use what we know about two and three dimensions to say a little more.

We saw that when we looked at the three-dimensional cube from the side, the blue center sphere was slightly bigger because of the extra dimension. If we lived in four-dimensional space, we could similarly look at the four-dimensional cube from the side so that it looked like a three-dimensional cube. From this angle, we would find that the corner spheres appear to overlap the center sphere. So the blue center sphere in four dimensions is even larger than the cube in three dimensions, which in turn is larger than the circle in two dimensions. However, the side lengths of the cubes are the same in each dimension!

Again, you could calculate the radius of the center sphere in the four-dimensional cube, and we would see that it is increasing. In fact, we can calculate the radius for dimension five, six and so on. An interesting thing happens when you get to dimension sixteen: The radius of the center sphere is exactly one half, which mean that it touches the sides of the (sixteen-dimensional) cube. In dimension seventeen, its radius is even bigger and it crashes through the sides of the cube. This is, of course, very counter-intuitive from the two- and three-dimensional pictures and it serves to illustrate the curse of dimensionality: While we can use our intuition from two and three dimensions to understand some aspects of higher dimensional geometry, there are also a lot of ways that our intuition can steer us wrong.

There is nothing we can do to eliminate this curse, but there are things we can do to minimize it. In particular, we should always look for ways to check that the methods that work well in low dimensions continue to be effective in higher dimensions. Having a solid understanding of the geometry behind different techniques can help us to both predict where these kinds of issues will arise, and find ways to double check that our results are still valid.

First, we need to look at regression from a slightly different perspective. When we originally fit a line to our data set, we treated the x and y coordinates interchangeably. We can also think of a line as a function that takes a value x and outputs a value y = cx + b for some pre-chosen parameters c and b that are determined by the regression algorithm. This makes it explicit that we are using the value of x to predict the value of y. Similarly, if we have a larger number of variables , we can describe a hyperplane by a function , where and are parameters that are calculated by the regression algorithm.

From this perspective, if we want to make regression more flexible, there seems to be a clear solution: We can replace linear function that defines a hyperplane with a more complicated function that has more parameters and thus can be better fit to the data. In two dimensions, for example, we could replace the line with the parabola . We now have three parameter (d, c and b) and the resulting graph will be allowed to curve up or down, as in the middle picture below.

We can again define a probability distribution concentrated along the parabola and use it to calculate the parameters that maximize the probability of a given data set, like we did with linear regression. But now we have three “dials” to tune with, instead of two, so we will be able to get a better value than with just a line. Adjusting each parameter moves the parabola in some way – shifting it up and down, left and right, or steeper and shallower, and the distribution will move with it. Rather than defining the distribution based on the distance from each data point to the curve, we’ll define it using the difference in the y-value between the data point and the point on the curve with the same x-value as the point. This turns out to be easier to work with when one actually writes the algorithm, and only changes the distribution slightly. (See the comments on the linear regression post for a discussion of the difference.) The standard approach is to take the Gaussian function of the difference squared, which gives us least squares regression, but there are other distributions one can use as well.

A parabola gives us more flexibility than a line, but it is still relatively rigid – It can only curve up or curve down. If, for example, our data follows an S-shaped curve, we will need at least a cubic curve () to describe it. As we try to model more complicated shapes, we need to add more terms and more parameters. We could also use functions besides polynomials. But this flexibility also means that we have to decide how much of the apparent complexity in the data set is the actual structure of the data and how much is noise. For example, in the picture on the right side of the Figure, the curve fits the data perfectly, but is much more complicated than is probably useful.

In particular, given any set of data points with distinct x-values, it is possible to find a polynomial function that exactly passes through each point, as in the picture on the right. This function will almost always be far more complicated than is needed to describe the data, and it’s an example of what’s called over-fitting: The distribution suggested by such a function does an excellent job of describing the existing data, but will do a lousy job of predicting new data. This is an issue that I touched on in the post on distributions. We will generally want to choose a polynomial somewhere in between – a function that has enough flexibility to capture the structure of the data, but restricted enough that it does not overfit.

But before we get into the problem of finding this perfect medium, we should consider general regression functions for higher dimensional data sets. If we have three variables x, y, z, we can replace the two-dimensional plane with a degree-two polynomial such as . This defines a two-dimensional shape called a paraboloid that curves in the and directions. This function has six parameters and we can “tune” them by defining a probability distribution based on the difference in the z-value between a data point the paraboloid. As always, we choose the parameters that maximize the probability of the data set with respect to the distribution defined by the parameters. If we increase the dimensions and/or increase the powers of the variables (such as ) then the number of parameters grows very quickly, but the model is always a codimension-one shape. (And again we can also use functions other than polynomials.) As the number of parameters increases, the model will fit the data points better, but it will also become more complicated and the risk of over fitting will increase.

So how do we select the correct number of parameters/level of complexity? There’s no hard and fast rule for this, and this is the area where data analysis take more the form of an arts than a science. If you have experience with a certain type of data (domain expertise) or have a rough theoretical model for the data, this can often guide your choice of a regression model. This is a problem where I suspect that a better understanding of the geometry could lead to some interesting new approaches. In general, though, the standard practice is to divide your data set into two sets, one called the training set T and the other called the evaluation set E. We will think of the training set as the existing data, and the evaluation set as the “new” data that we want the regression model to predict. In the pictures below, we’ve left the training set colored blue, but changed the evaluation set to orange.

If we run regression on the training set with a relatively simple curve, as on the left, the curve still stays fairly close to the evaluation set. However, if we run regression with the more flexible curve, but using only the data points from training set then the resulting curve (on the right) passes through every training point, but is farther from the evaluation points.

In general, for each choice of a type of function (line/parabola/etc.), we’ll choose the parameters by maximizing the probability of the training set T for the resulting distribution. This may be overfit or underfit, but we can check this by calculating the probability of the evaluation set E with respect to each resulting distribution. The one with the best value for E is the least likely to be overfit. So even if we do this for a higher dimensional data set in which we can’t actually see the resulting “curves”, the scores for the evaluation set would suggest which curves suffer from over-fitting. If we do this for many different types of functions and choose the one that does the best job of predicting the evaluation data, we can be relatively confident that this model is a good median between overfitting and underfitting.

In two-dimensional linear regression, the general form for a model is a distribution concentrated along a line. A line is determined by two parameters – its slope and it y-intercept – and we want to find the parameters that determine the best fit line for a given set of points. We know that the data points probably won’t all fall right on any one line, so there will always be some error. This is why we have to allow for a blurry probability distribution, as in the last post.

For any given line, we can define a distribution that is equal to one along the line and decreases as we move away from the line. In particular, the probability will be defined by the Gaussian function ( where is the distance), so that as we move away from the line, the probability will follow a bell curve shown in the Figure above. The right side shows a graph of the distribution function – The height of the graph is the probability at each point. This looks like a bell curve extended along the line.

The regression algorithm picks out the one line, out of all possible lines, for which the data points best fit the corresponding distribution. But how do we determine how good the fit is? A given distribution assigns to each point a probability between 0 and 1, indicated below by how dark the point is. You can think of this as the probability that the given point will be randomly chosen. The probability that a collection of points would be chosen at random is the product of their individual probabilities. The regression algorithm chooses the distribution that assigns the highest probability to the given data points. In other words, it chooses the distribution in which the data points are overall in the darkest possible parts of the distribution.

For example, in the figure above, the regression algorithm would choose a line something like the one on the left rather than the one on the right for the blue data points. You can think of this as adjusting two parameters – the slope of the line and its -intercept – until the probability value is maximized. This is kind of like tuning an old-fashioned analog radio: As you move the knob back and forth, the signal gets stronger and weaker and you stop when the signal is as strong as possible. You might also imagine printing the distribution onto a sheet of clear plastic, placing it over the data points and moving it around (translating and rotating) until you find the position where the points are in the darkest region possible.

In practice, the regression algorithm doesn’t actually try all the different parameter values because we can switch to a log scale (which turns multiplying the probabilities into adding the sum of the squared distances) and then apply some tricks from calculus to directly calculate the ideal line. If you’re interested in the details of this, there are plenty of good statistics books out there. Since I want to focus on the geometric intuition, I’m going to skip over it here.

There are two ways to generalize two-dimensional linear regression: We can change two-dimensional to higher dimensional and/or we can change the line to a more flexible shape. I’ll introduce more flexible shapes in the next post, so for now lets consider higher dimensional linear regression.

The usual goal of regression is to predict the value of one variable/dimension based on the other value(s) of the other variable(s). If the two dimensions represent height and age of a collection of trees, then the best fit line will allow you estimate the age of a new tree based on its height. This works because on the line, each x-value determines a unique y-value. In three dimensions, a two-dimensional plane has a similar property: If we know the values of any two of x, y or z then we can predict the third. Three-dimensional regression therefore involves fitting a plane to a data set with a distribution that is equal to one on the plane and follows a bell curve as we get farther from the plane. For example, if the data set records the height, width and age of a collection of trees, the best fit plane will allow you to predict the age of a tree based on its height and width (or predict the height of a tree based on its width and age, etc.)

Notice that for both two-dimensional and three-dimensional linear regression, the distribution is concentrated near a shape whose dimension is one less than the whole space (a one-dimensional line or a two-dimensional plane, respectively). The technical term for the dimension of a space minus the dimension of a shape in that space is the co-dimension of the shape. So in two- and three-dimensional linear regression, we’re looking for a shape whose codimension is equal to one. In general, in a codimension-one shape, the value of one variable will be determined by the other variables.

In four-dimensional space, a shape with co-dimension one will be three-dimensional. So four-dimensional linear regression (i.e. regression with four variables) involves distributions concentrated near three-dimensional hyperplanes. There’s essentially no hope of visualizing this, but the two- and three-dimensional cases should give you some idea of how it works. Each point is defined by three variables, say x, y, z and w and once you’ve found the best fit hyperplane, if you know any three of the values then you can predict the fourth. Higher dimensional linear regression is similar, and in each dimension the goal is to find the best distribution concentrated near a hyperplane whose dimension is one less than that of the whole space.

References: In the comments, Saeid suggested that I add references. I don’t know the Statistics literature very well, but the following books look helpful, based on Amazon.com. If there are any books you like for their description of regression, please let me know in the comments.

• The Geometry of Multivariate Statistics – Thomas D. Wickens
• Statistical Methods: A Geometric Primer – David Saville & Graham Wood
• Statistical Methods: The Geometric Approach – David Saville & Graham Wood
• Advanced Data Analysis from an Elementary Point of View – Cosma Rohilla Shalizi (with a link to a PDF of the first draft)
• The Elements of Statistical Learning – Trevor Hastie, Robert Tibshirani and Jerome Friedman