Flash and JavaScript are required for this feature.
Download the video from iTunes U or the Internet Archive.
Topics covered: Validating simulation results, curve fitting, linear regression
Instructor: Prof. Eric Grimson, Prof. John Guttag
21: Validating Simulation R...
OPERATOR: The following content is provided under a Creative Commons license. Your support will help MIT OpenCourseWare continue to offer high quality educational resources for free. To make a donation or view additional materials from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu.
PROFESSOR: So let's start. I have written a number on the board here. Anyone want to speculate what that number represents? Well, you may recall at the end of the last lecture, we were simulating pi, and I started up running it with a billion darts. And when it finally terminated, this was the estimate of pi it gave me with a billion. Not bad, not quite perfect, but still pretty good. In fact when I later ran it with 10 billion darts, which took a rather long time to run, didn't do much better. So it's converging very slowly now near the end.
When we use an algorithm like that one to perform a Monte Carlo simulation, we're trusting, as I said, that fate will give us an unbiased sample, a sample that would be representative of true random throws. And, indeed in this case, that's a pretty good assumption. The random number generator is not truly random, it's what's called pseudo-random, in that if you start it with the same initial conditions, it will give you the same results. But it's close enough for, at least for government work, and other useful projects. We do have to think about the question, how many samples should we run? Was a billion darts enough? Now since we sort of all started knowing what pi was, we could look at it and say, yeah, pretty good. But suppose we had no clue about the actual value of pi. We still have to think about the question of how many samples? And also, how accurate do we believe our result is, given the number of samples? As you might guess, these two questions are closely related. That, if we know in advance how much accuracy we want, we can sometimes use that to calculate how many samples we need.
But there's still always the issue. It's never possible to achieve perfect accuracy through sampling. Unless you sample the entire population. No matter how many samples you take, you can never be sure that the sample set is typical until you've checked every last element. So if I went around MIT and sampled 100 students to try and, for example, guess the fraction of students at MIT who are of Chinese descent. Maybe 100 students would be enough, but maybe I would get unlucky and draw the wrong 100. In the sense of, by accident, 100 Chinese descent, or 100 non-Chinese descent, which would give me the wrong answer. And there would be no way I could be sure that I had not drawn a biased sample, unless I really did have the whole population to look at.
So we can never know that our estimate is correct. Now maybe I took a billion darts, and for some reason got really unlucky and they all ended up inside or outside the circle. But what we can know, is how likely it is that our answer is correct, given the assumptions. And that's the topic we'll spend the next few lectures on, at least one of the topics. It's saying, how can we know how likely it is that our answer is good. But it's always given some set of assumptions, and we have to worry a lot about those assumptions. Now in the case of our pi example, our assumption was that the random number generator was indeed giving us random numbers in the interval 0 to 1. So that was our underlying assumption. Then using that, we looked at a plot, and we saw that after time the answer wasn't changing very much. And we use that to say, OK, it looks like we're actually converging on an answer. And then I ran it again, with another trial, and it converged again at the same place. And the fact that that happened several times led me to at least have some reason to believe that I was actually finding a good approximation of pi.
That's a good thing to do. It's a necessary thing to do. But it is not sufficient. Because errors can creep into many places. So that kind of technique, and in fact, almost all statistical techniques, are good at establishing, in some sense, the reproduce-ability of the result, and that it is statistically valid, and that there's no error, for example, in the way I'm generating the numbers. Or I didn't get very unlucky. However, they're other places other than bad luck where errors can creep in.
So let's look at an example here. I've taken the algorithm we looked at last time for finding pi, and I've made a change. You'll remember that we were before using 4 as our multiplier, and here what I've done is, just gone in and replaced 4 by 2. Assuming that I made a programming error. Now let's see what happens when we run it. Well, a bad thing has happened. Sure enough, we ran it and it converged, started to converge, and if I ran 100 trials each one would converge at roughly the same place. Any statistical test I would do, would say that my statistics are sound, I've chosen enough samples, and for some accuracy, it's converting. Everything is perfect, except for what? It's the wrong answer. The moral here, is that just because an answer is statistically valid, does not mean it's the right answer.
And that's really important to understand, because you see this, and we'll see more examples later, not today, but after Thanksgiving, comes up all the time in the newspapers, in scientific articles, where people do a million tests, do all the statistics right, say here's the answer, and it turns out to be completely wrong. And that's because it was some underlying assumption that went into the decision, that was not true. So here, the assumption is, that I've done my algebra right for computing pi based upon where the darts land. And it turns out, if I put 2 here, my algebra is wrong. Now how could I discover this? Since I've already told you no statistical test is going to help me. What's the obvious thing I should be doing when I get this answer? Somebody? Yeah?
STUDENT: [INAUDIBLE]
PROFESSOR: Exactly. Checking against reality. I started with the notion that pi had some relation to the area of a circle. So I could use this value of pi, draw a circle with a radius. Do my best to measure the area. I wouldn't need to get a very good, accurate measurement, and I would say, whoa, this isn't even close. And that would tell me I have a problem. So the moral here is, to check results against physical reality. So for example, the current problem set, you're doing a simulation about what happens to viruses when drugs are applied. If you were doing this for a pharmaceutical company, in addition to the simulation, you'd want to run some real experiments. And make sure that things matched. OK, what this suggests, is that we often use simulation, and other computational techniques, to try and model the real world, or the physical world, in which we all live. And we can use data to do that. I now want to go through another set of examples, and we're going to look at the interplay of three things: what happens when you have data, say from measurements, and models that at least claim to explain the data. And then, consequences that follow from the models. This is often the way science works, its the way engineering works, we have some measurements, we have a theory that explains the measurements, and then we write software to explore the consequences of that theory. Including, is it plausible that it's really true?
So I want to start, as an example, with a classic chosen from 8.01. So I presume, everyone here has taken 8.01? Or in 8.01? Anyone here who's not had an experience with 801? All right, well. I hope you know about springs, because we're going to talk about springs. So if you think about it, I'm now just talking not about springs that have water in them, but springs that you compress, you know, and expand, and things like that. And there's typically something called the spring constant that tells us how stiff the spring is, how much energy it takes to compress this spring. Or equivalently, how much pop the spring has when you're no longer holding it down. Some springs are easy to stretch, they have a small spring constant. Some strings, for example, the ones that hold up an automobile, suspension, are much harder to stretch and compress.
There's a theory about them called Hooke's Law. And it's quite simple. Force, the amount of force exerted by a spring, is equal to minus some constant times the distance you have compressed the spring. It's minus, because the force is exerted in an opposite direction, trying to spring up. So for example, we could look at it this way. We've got a spring, excuse my art here. And we put some weight on the spring, which has therefore compressed it a little bit. And the spring is exerting some upward force. And the amount of force it's exerting is proportional to the distance x. So, if we believe Hooke's Law, and I give you a spring, how can we find out what this constant is? Well, we can do it by putting a weight on top of the spring. It will compress the spring a certain amount, and then the spring will stop moving. Now gravity would normally have had this weight go all the way down to the bottom, if there was no spring. So clearly the spring is exerting some force in the upward direction, to keep that mass from going down to the table, right? So we know what that force is there. If we compress the spring to a bunch of different distances, by putting, say, different size weights on it, we can then solve for the spring constant, just the way, before, we solved for pi.
So it just so happens, not quite by accident, that I've got some data from a spring. So let's look at it. So here's some data taken from measuring a spring. This is distance and force, force computed from the mass, basically, right? Because we know that these have to be in balance. And I'm not going to ask you to in your head estimate the constant from these, but what you'll see is, the format is, there's a distance, and then a colon, and then the force. Yeah?
STUDENT: [INAUDIBLE]
PROFESSOR: Ok, right, yes, thank you. All right, want to repeat that more loudly for everyone?
STUDENT: [INAUDIBLE]
PROFESSOR: Right, right, because the x in the equation -- right, here we're getting an equilibrium. OK, so let's look at what happens when we try and examine this. We'll look at spring dot pi. So it's pretty simple. First thing is, I've got a function that reads in the data and parses it. You've all done more complicated parsing of data files than this. So I won't belabor the details. I called it get data rather than get spring data, because I'm going to use the same thing for a lot of other kinds of data. And the only thing I want you to notice, is that it's returning a pair of arrays. OK, not lists. The usual thing is, I'm building them up using lists, because lists have append and arrays don't, and then I'm converting them to arrays so I can do matrix kinds of operations on them.
So I'll get the distances and the forces. And then I'm just going to plot them, and we'll see what they look like. So let's do that. There they are. Now, if you believe Hooke's Law, you could look at this data, and maybe you wouldn't like it. Because Hooke's Law implies that, in fact, these points should lie in a straight line, right? If I just plug in values here, what am I going to get? A straight line, right? I'm just multiplying k times x. But I don't have a straight line, I have a little scatter of points, kind of it looks like a straight line, but it's not. And why do you think that's true? What's going on here? What could cause this line not to be straight? Have any you ever done a physics experiment? And when you did it, did your results actually match the theory that your high school teacher, say, explained to you. No, and why not. Yeah, you have various kinds of experimental or measurement error, right? Because, when you're doing these experiments, at least I'm not perfect, and I suspect at least most of you are not perfect, you get mistakes. A little bit of error creeps in inevitably. And so, when we acquired this data, sure enough there was measurement error. And so the points are scattered around.
This is something to be expected. Real data almost never matches the theory precisely. Because there usually is some sort of experimental error that creeps into things. So what should we do about that? Well, what usually people do, when they think about this, is they would look at this data and say, well, let me fit a line to this. Somehow, say, what would be the line that best approximates these points? And then the slope of that line would give me the spring constant.
So that raises the next question, what do I mean by finding a line that best fits these points? How do we, fit, in this case, a line, to the data? First of all, I should ask the question, why did I say let's fit a line? Maybe I should have said, let's fit a parabola, or let's fit a circle? Why should I had said let's fit a line. Yeah?
STUDENT: [INAUDIBLE]
PROFESSOR: Well, how do I know that the plot is a linear function? Pardon? Well, so, two things. One is, I had a theory. You know, I had up there a model, and my model suggested that I expected it to be linear. And so if I'm testing my model, I should and fit a line, my theory, if you will. But also when I look at it, it looks kind of like a line. So you know, if I looked at it, and it didn't look like a line, I might have said, well, my model must be badly broken. So let's try and see if we can fit it. Whenever we try and fit something, we need some sort of an objective function that captures the goodness of a fit. I'm trying to find, this is an optimization problem of the sort that we've looked at before. I'm trying to find a line that optimizes some objective function. So a very simple objective function here, is called the least squares fit. I want to find the line that minimizes the sum of observation sub i, the i'th data point I have, minus what the line, the model, predicts that point should have been, and then I'll square it. So I want to minimize this value. I want to find the line that gives me the smallest value for this.
Why do you think I'm squaring the difference? What would happen if I didn't square the difference? Yeah? Positive and negative errors might cancel each other out. And in judging the quality of the fit, I don't really care deeply -- you're going to get very fat the way you're collecting candy here -- I don't care deeply whether the error is, which side, it is, just that it's wrong. And so by squaring it, it's kind of like taking the absolute value of the error, among other things. All right, so if we look at our example here, what would this be? I want to minimize, want to find a line that minimizes it. So how do I do that? I could easily do it using successive approximation, right? I could choose a line, basically what I am, is I'm choosing a slope, here, right? And, I could, just like Newton Raphson, do successive approximation for awhile, and get the best fit. That's one way to do the optimization. It turns out that for this particular optimization there's something more efficient. You can actually, there is a closed form way of attacking this, and I could explain that, but in fact, I'll explain something even better. It's built into Pylab.
So Pylab has a function built-in called polyfit. Which, given a set of points, finds the polynomial that gives you the best least squares approximation to those points. It's called polynomial because it isn't necessarily going to be first order, that is to say, a line. It can find polynomials of arbitrary degree. So let's look at the example here, we'll see how it works. So let me uncomment it. So I'm going to get k and b equals Pylab dot polyfit here. What it's going to do is, think about a polynomial. I give you a polynomial of degree one, you have all learned that it's a x plus b, b is the constant, and x is the single variable. And so I multiply a by x and I add b to it, and as I vary x I get new values. And so polyfit, in this case, will take the set of points defined by these two arrays and return me a value for a and a value for b. Now here I've assigned a to k, but don't worry about that. And then, I'm gonna now generate the predictions that I would get from this k and b, and plot those.
So let's look at it. So here it said the k is 31.475, etc., and it's plotted the line that it's found. Or I've plotted the line. You'll note, a lot of the points don't lie on the line, in fact, most of the points don't lie on the line. But it's asserting that this is the best it can do with the line. And there's some points, for example, up here, that are kind of outliers, that are pretty far from the line. But it has minimized the error, if you will, for all of the points it has. That's quite different from, say, finding the line that touches the most points, right? It's minimizing the sum of the errors. Now, given that I was just looking for a constant to start with, why did I bother even plotting the data? I happen to have known before I did this that polyfit existed, and what I was really looking for was this line. So maybe I should have just done the polyfit and said here's k and I'm done. Would that have been a good idea? Yeah?
STUDENT: You can't know without seeing the actual data how well it's actually fitting it.
PROFESSOR: Right. Exactly right. That says, well how would I know that it was fitting it badly or well, and in fact, how would I know that my notion of the model is sound, or that my experiment isn't completely broken? So always, I think, always look at the real data. Don't just, I've seen too many papers where people show me the curve that fits the data, and don't show me the data, and it always makes me very nervous. So always look at the data, as well as however you're choosing to fit it.
As an example of that, let's look at another set of inputs. This is not a spring. It's the same get data function as before, ignore that thing at the top. I'm going to analyze it and we'll look at it. So here I'm plotting the speed of something over time. So I plotted it, and I've done a least squares fit using polyfit just as before to get a line, and I put the line vs. the data, and here I'm a little suspicious. Right, I fit a line, but when I look at it, I don't think it's a real good fit for the data. Somehow modeling this data as a line is probably not right. A linear model is not good for this data. This data is derived from something, a more complex process. So take a look at it, and tell me what order were calling of polynomial do you think might fit this data? What shape does this look like to you? Pardon?
STUDENT: Quadratic.
PROFESSOR: Quadratic, because the shape is a what? It's a parabola. Well, I don't know if I dare try this one all the way to the back. Ooh, at least I didn't hurt anybody. All right, fortunately it's just as easy to fit a ravel parabola as a line. So let's look down here. I've done the same thing, but instead of passing it one, as I did up here as the argument, I'm passing it two. Saying, instead of fitting a polynomial of degree one, fit a polynomial of degree two. And now let's see what it looks like. Well, my eyes tell me this is a much better fit than the line. So again, that's why I wanted to see the scatter plot, so that I could at least look at it with my eyes, and say, yeah, this looks like a better fit. All right, any question about what's going on here?
What we've been looking at is something called linear regression. It's called linear because the relationship of the dependent variable y to the independent variables is assumed to be a linear function of the parameters. It's not because it has to be a linear function of the value of x, OK? Because as you can see, we're not getting a line, we're getting a parabola. Don't worry about the details, the point I want to make is, people sometimes see the word linear regression and think it can only be used to find lines. It's not so. So when, for example, we did the quadratic, what we had is y equals a x squared plus b x plus c. The graph vs. x will not be a straight line, right, because I'm squaring x. But it is, just about, in this case, the single variable x. Now, when I looked at this, I said, all right, it's clear that the yellow curve is a better fit than the red. It's a red line. But that was a pretty informal statement.
I can actually look at this much more formally. And we're going to look at something that's the statisticians call r squared. Which in the case of a linear regression is the coefficient of determination. Now, this is a big fancy word for something that's actually pretty simple. So what r squared its going to be, and this is on your handout, is 1 minus e e over d v. So e e is going to be the errors in the estimation. So I've got some estimated values, some predicted values, if you will, given to me by the model, either the line or the parabola in this case. And I've got some real values, corresponding to each of those points, and I can look at the difference between the 2 And that will tell me how much difference there is between the estimated data and the, well, between the predicted data and the measured data, in this case. And then I want to divide that by the variance in the measured data. The data variance. How broadly scattered the measured points are. And I'll do that by comparing the mean of the measured data, to the measured data. So I get the average value of the measured data, and I look at how different the points I measure are.
So I just want to give to you, informally, because I really don't care if you understand all the math. What I do want you to understand, when someone tells you, here's the r squared value, is, informally what it really is saying. It's attempting to capture the proportion of the response variation explained by the variables in the model. In this case, x. So you'll have some amount of variation that is explained by changing the values of the variables. So if, actually, I'm going to give an example and then come back to it more informally. So if, for example, r squared were to equal 0.9, that would mean that approximately 90 percent of the variation in the variables can be explained by the model. OK, so we have some amount of variation in the measured data, and if r squared is 0.9, it says that 90 percent can be explained by the models, and the other 10 percent cannot. Now, that other 10 percent could be experimental error, or it could be that, in fact, you need more variables in the model. That there are what are called lurking variables. I love this term. A lurking variable is something that actually effects the result, but is not reflected in the model.
As we'll see a little bit later, this is a very important thing to worry about, when you're looking at experimental data and you're building models. So we see this, for example, in the medical literature, that they will do some experiment, and they'll say that this drug explains x, or has this affect. And the variables they are looking at are, say, the disease the patient has, and the age of the patient. Well, maybe the gender of the patient is also important, but it doesn't happen to be in the model. Now, if when they did a fit, it came out with 0.9, that says at worst case, the variables we didn't consider could cause a 10 percent error. But, that could be big, that could matter a lot. And so as you get farther from 1, you ought to get very worried about whether you actually have all the right variables. Now you might have the right variables, and just experiment was not conducted well, But it's usually the case that the problem is not that, but that there are lurking variables. And we'll see examples of that.
So, easier to read than the math, at least by me, easier to read than the math, is the implementation of r square. So it's measured and estimated values, I get the diffs, the differences, between the estimated and the measured. These are both arrays, so I subtract 1 array from the other, and then I square it. Remember, this'll do an element-wise subtraction, and then square each element. Then I can get the mean, by dividing the sum of the array measured by the length of it. I can get the variance, which is the measured mean minus the measured value, again squared. And then I'll return 1 minus this. All right? So, just to make sure we sort of understand the code, and the theory here as well, what would we get if we had absolutely perfect prediction? So if every measured point actually fit on the curb predicted by our model, what would r square return? So in this case, measured and estimated would be identical. What gets return by this? Yeah, 1. Exactly right. Because when I compute it, it will turn out that these two numbers will be the, I'll get 0, 1 minus 0 is 0, right? Because the differences will be zero. OK?
So I can use this, now, to actually get a notion of how good my fit is. So let's look at speed dot pi again here, and now I'm going to uncomment these two things, where I'm going to, after I compute the fit, I'm going to then measure it. And you'll see here that the r squared error for the linear fit is 0.896, and for the quadratic fit is 0.973. So indeed, we get a much better fit here. So not only does our eye tell us we have a better fit, our more formal statistical measure tells us we have a better fit, and it tells us how good it is. It's not a perfect fit, but it's a pretty good fit, for sure.
Now, interestingly enough, it isn't surprising that the quadratic fit is better than the linear fit. In fact, the mathematics of this should tell us it can never be worse. How do I know it can never be worse? That's just, never is a really strong word. How do I know that? Because, when I do the quadratic fit, if I had perfectly linear data, then this coefficient, whoops, not that coefficient, wrong, this coefficient, could be 0. So if I ask it to do a quadratic fit to linear data, and the a is truly perfectly linear, this coefficient will be 0, and my model will turn out to be the same as the linear model. So I will always get at least as good a fit.
Now, does this mean that it's always better to use a higher order polynomial? The answer is no, and let's look at why. So here what I've done is, I've taken seven points, and I've generated, if you look at this line here, the y-values, for x in x vals, points dot append x plus some random number. So basically I've got something linear in x, but I'm perturbing, if you will, my data by some random value. Something between 0 and 1 is getting added to things. And I'm doing this so my points won't lie on a perfectly straight line. And then we'll try and fit a line to it. And also, just for fun, we'll try and fit a fifth order polynomial to it. And let's see what we get. Well, there's my line, and there's my fifth order polynomial. Neither is quite perfect, but which do you think looks like a closer fit? With your eye. Well, I would say the red line, the red curve, if you will, is a better fit, and sure enough if we look at the statistics, we'll see it's 0.99, as opposed to 0.978. So it's clearly a closer fit.
But that raises the very important question: does closer equal better, or tighter, which is another word for closer? And the answer is no. It's a tighter fit, but it's not necessarily better, in the sense of more useful. Because one of the things I want to do when I build a model like this, is have something with predictive power. I don't really necessarily need a model to tell me where the points I've measured lie, because I have them. The whole purpose of the model is to give me some way to predict where unmeasured points would lie, where future points would lie. OK, I understand how the spring works, and I can guess where it would be if things I haven't had the time to measure, or the ability to measure. So let's look at that. Let's see, where'd that figure go. It's lurking somewhere. All right, we'll just kill this for now. So let's generate some more points, and I'm going to use exactly the same algorithm. But I'm going to generate twice as many points. But I'm only fitting it to the first half. So if I run this one, figure one is what we looked at before. The red line is fitting them a little better. But here's figure two. What happens when I extrapolate the curve to the new points? Well, you can see, it's a terrible fit. And you would expect that, because my data was basically linear, and I fit in non-linear curve to it. And if you look at it you can see that, OK, look at this, to get from here to here, it thought I had to take off pretty sharply. And so sure enough, as I get new points, the prediction will postulate that it's still going up, much more steeply than it really does. So you can see it's a terrible prediction.
And that's because what I've done is, I over-fit the data. I've taken a very high degree polynomial, which has given me a good close fit, and I can always get a fit, by the way. If I choose a high enough degree polynomial, I can fit lots and lots of data sets. But I have reason to be very suspicious. The fact that I took a fifth order polynomial to get six points should make me very nervous. And it's a very important moral. Beware of over-fitting. If you have a very complex model, there's a good chance that it's over-fit. The larger moral is, beware of statistics without any theory. You're just cranking away, you get a great r squared, you say it's a beautiful fit. But there was no real theory there. You can always find a fit. As Disraeli is alleged to have said, there are three kinds of lies: lies, damned lies, and statistics. And we'll spend some more time when we get back from Thanksgiving looking at how to lie with statistics. Have a great holiday, everybody.
Free Downloads
Video
- iTunes U (MP4 - 118MB)
- Internet Archive (MP4 - 118MB)
Free Streaming
Subtitle
- English - US (SRT)