Chapter 21.1 - Data Assimilation

0:08 >> We begin a new chapter on what's called data simulation, which is a really important technique for us when we start thinking about building models and collecting data and using models and data jointly. In fact, data simulation is really why we've made such big advancements in weather prediction and making longer time forecasts because what we're doing there is using our data and our models jointly together to make the best predictions possible. So, it's exactly what we're going to do here, and I want to show you the basic framework for it because it's essentially something that is specially as you do something in practice where you can collect real data, the assimilation procedure is going to be really important to taking the model, which is an approximation to your physics, and improving on it.

0:55 So, really what we're going to do is thinking about using models and data jointly, and data simulation is in a framework, a mathematical framework, to help you do that. So, let's first of all start with the model. Typically, when you're taking a system, you have some kind of physics-based model and generically represented here by a dynamical system. d y d t equals to f, so f is the prescription for our governing equations. This is a generic framework, right? So, this could be something like weather models where I'm predicting global circulation, it could be electromagnetism, it could be fluid dynamics. Whatever it is, the idea is that I have a model for making a future state prediction.

1:37 And in that model, the way I run it is I give you an initial condition. So, the model is not just prescribing this f that we typically build from physics-type models, but also prescribing the initial condition. So, this is how we've actually engineered a a of science and engineering practice in in in the past as we've always taken a model like this, done simulations, and started to learn how to build new things or how to make predictions in an accurate and meaningful way.

2:11 But, both your model is wrong and your initial conditions are often only measured to a certain degree of accuracy. So, I want to talk about these two quantities here, Q1 and Q2. So, what this is is to say, "Look, whatever model you've picked, there's a couple things that can happen. Either you've potentially neglected some of the physics, so what you have here in the F is sort of, let's say, the dominant physics of the system, but when you actually, you know, it's an approximation to what you think the physics is. Q1 represents everything else that's there, including potentially noise-driven stochastic terms underneath at small scales or bigger scales that you haven't correctly modeled with the function F.

2:54 Okay? So, Q1 is a generic term basically telling you your model and reality mismatch. Your model is not the truth. Your model is an approximation to the truth. Q1 is, in some sense, you're not going to have access to Q1. It's just an acknowledgement, in some sense, that there are physics pieces missing to your model. Q2 is essentially your inability to exactly measure the initial state.

3:25 Why not? So, for instance, you're going to take, right, you're going to start to simulate something. You have to give it an initial condition. For instance, in weather prediction, you have to say, "From the what the weather is right now, what does the future look like?" So, you would have to have an initial condition, but your measurements only going to approximate this Y0, and Q2 tells you the mismatch between your approximation and reality. So, really, this is getting to now beyond sort of building quantitative models for reality to now trying to tie those models directly to reality and acknowledging the fact that all of our models are inadequate for monitoring modeling true reality. Okay?

4:09 So, that's on the model side. Let's talk about on the data side. The data side is has a similar issue which is typically measurements, you know, you might take you can always write a measurement here as some function G is equal to zero. Could for instance be the pressure. So, you could say the pressure is 10. So, you'd say pressure minus 10 is equal to zero. So, that's what the G would be. But Q3 represents for us the noise in that measurement. In other words, your measuring device itself has a certain accuracy and the minute that you lose some of that accuracy or that you have noise in those sensors, Q3 is aimed at basically identifying the imperfections or the measurement error or noise in your sensor. This is a common thing. It There's no way you can make really a perfect sensor. And so, your sensor has when you take a measurement, it has a close value but not a perfect value. And what we want to think about with Q3 is what's the error distribution from my noisy measurements. And it could simply be that your sensor has noisy measurements and so that would be Q3.

5:19 And that's typically what we're going to assume in all of this architecture. So, overall, look at the problem that we have when we have both model and data jointly. In a lot of this book, we've either focused on I have models and I simulate them or we've built models directly from data. And the data simulation framework acknowledges two things. One, that your model, which includes your you know, missing physics or discrepancies that you might have with real physics as well as your inability to prescribe initial conditions this Q1 and Q2 acknowledges that those are imperfect in that if you have these errors in your model or your in your initial condition there can be a problem.

6:02 And in here it also acknowledges the fact the second fact that whatever measurement you take of the system is not perfect. In fact, it is inaccurate with some error. And so what we want to do is try to think about how do the start using these together to make the best predictions. We're going to focus here on making for instance forecasts of a dynamical system where we take our model to make the forecast but we take the measurement to pin that model more accurately to reality. That's the whole idea of data simulation.

6:35 And in doing so, we define this introduces quadratic form. In some sense what we're going to do here, check out what we've got here. We've got the Q1 part of the form here. So this is integrated over the trajectory of the dynamics, okay? So Q1, remember, is right here. This is my error in my dynamical system. So it's a function of T. So as I solve this system from zero to some capital T or some time range, Q1 accumulates over that.

7:06 Q2 is just the initial condition. Q3 is at the measurement. So the way I'm going to think about this is Q1, I integrate it as a quadratic form with some weighting matrix on my W1 here where now I integrate over the trajectory that I would have from time zero to T, okay? That is going to be essentially thinking about the error of my trajectory. This is the error or the variance because it's a quadratic form of essentially Q2 which is my variance in my error in my condition, and this is the variance in my sensor measurements. So, all together, we can think about these different pieces all contributing to my overall error in making predictions. And so, what we're going to define is this quadratic form because what we're going to do is try to find a way to minimize the variance in the error and use the measurements and the errors together to get the best predictions possible.

8:09 So, that's the interpretation of what this variance means. So, variance of the trajectories, the variance of the initial conditions, the variance of the measurements. So, we're going to use this, and we're going to try to optimize over this in order to reduce the error as much as we can. Okay. This leads us directly to what's called Bayes' theorem. And what we're going to think about this is generically Bayes' theorem is stated as a probability, usually a conditional probability, the probability of the variable X, what's its probability given that I have Y? In other words, if I'm not X isn't just I just don't look at the probability of X alone, I give it the probability of this value of X given that the value of Y is something in particular. This here can be written as the probability of Y given X times the probability of X divided by probability of Y. Now, the probability of Y itself is just a normalization there, so we're actually going to ignore it. And what we really want to focus is on these two terms here. And essentially, what we're going to think about here, X is going to be our model distribution, and okay, that's the P of X. So, we can look at when I do, for instance, when I run my model, what I can look at is especially under perturbations or under different initial conditions, I can ask the question, if I run those, what's the distribution of my model?

9:36 And then here, this is now is the distribution Y is my observation given my model X. So, we can compute these two because we can actually make these observation versus here, this is harder for us to do which is the probability of the model given the observation. We instead just this here we can directly observe this distribution and we can compute this from our model. We're going to ignore this cuz it's just a normalization.

10:06 And so, typically what's done in the simplest instances is you basically assume that these distributions are Gaussian. unless you have better information, that's typically what's always done in these statistical models is you say, "Well, I've got this model. What how is the distribution?" The first thing to try always is a Gaussian distribution. And so, what we assume here for instance, here's the distribution of the model around right center around some X naught which is my initial condition as it were. And here it's a distribution there with some variant sigma naught. And then this here is the distribution of my observation centered around X itself. So, this is Y is my observation, X is my model, and there's some variance Y there. And so, I'm going to assume that those are the distributions that I have working with.

10:56 So, the conditional probability, right, is this P YX. I just put those Gaussians together and this is what it looks like right here. Okay? Now, the nice thing is we're going to put this distribution into our quadratic form that we have. In the quadratic form, we're going to take the log of this probability distribution. In fact, I'm going to take the negative log of this. And I'm going to add this log of C3 cuz when I take the log of this, the log of the product of C3 constant that's in front which is related to the P of Y, it's just a constant I don't care about.

11:33 So, but if I take the log of this, I can take the log of I can break the log up into an addition of logs. And the log C3 is going to cancel with this. And the only reason I do that, I get rid of a normalization term. And then the log of the rest of this, basically, when I do that, these are now sums of two exponentials. The log of the exponential just brings out, in fact, just this quadratic right here. And in fact, oftentimes when we think about these quadratic forms in the log, we put the log in there explicitly cuz we are going to assume we're going to working with Gaussian distribution, and the exponent goes away directly out of that.

12:11 So, this is the way we construct our variance model, or essentially our quadratic form here, based upon these. And the nice thing is, we've constructed it, and we have this very simple result about what we have here. So, let's draw a picture of what we've actually done. So, what we're trying to do is make a best guess of what my solution should look like. So, I'm going to give you three things here. So, here's this bump here. This is let's say the P of X. This is my model distribution. I know my model's not right, but here's its distribution.

12:49 P of Y of X, this is my observational distribution given these X's that I have. In other words, when I make these predictions of the model, I take an actual actual observation of the model, and that's its distribution. Let's say this Gaussian. So, how should I use these together? What Bayes theorem says is, well, I can basically compute this here, which I've shown you, which is what's the probability of an X given measurement Y, is the same as the probability of Y given X, probability of X. And so, I take the product of these two Gaussians, normalized, and that's the bolded line.

13:26 In other words, my prediction is essentially using the model and the observations jointly because if I if if my observations are here, it's awkward to pick a value over here, even though my model said it could be here, when in fact, my observations say it can't be there. Okay? So, I'm biasing towards using these measurements, but I also want to use the model itself to try to help make predictions. And so, when we look at the product of these, I get this very narrow Gaussian here, which is essentially using Bayes' theorem.

14:05 And in fact, what we could do then is say, well, where is in fact the most likely value of this to be? And what I can do is just compute the derivative of that variational function with respect to X, set it to zero. And when you do that, here's what you get. This is my solution. What you would actually do if you only had the model is you would say, X is just run your model, and you could just take the the mean the you know the the maximum here of that Gaussian, that's your prediction.

14:37 Or you could just use your data, which would be Y. But here, what you've done is use use the two together, and this is in fact your prediction that weights the model and the observations together. Okay? In order to find yourself a solution here where the model itself is informed by the observation. So, that is the prediction value, X-bar. Here's the observation, here is your prediction from the model, here are these variance terms, which tell you essentially how fat these Gaussians are and how they come together into this Gaussian here.

15:14 We can also compute, just do a little bit of linear algebra around this, and compute the variance. So, the variance, sigma bar squared, of this one here is given by either here, you can either write in terms of sigma naught or sigma y, the variance of the measurement or the variance of the model, and K, either one will work, but the most important thing to note is that both of these are less than either the variance of your model, sigma naught squared, or the variance of your measurement, sigma y squared. So, when you do this assimilation trick, not only are you getting an improved prediction, but your variance is shrinking.

15:53 Okay? From either your measurement or from your model. So, that's great news. You're kind of getting a double win here, a better prediction, lower and you know, lowering the uncertainty, all from just taking this here and making these predictions. So, this is the way that assimilation works, and you can see it's a very powerful framework in a system where you both have a model and you have access to making measurements. And it should essentially always be used if you can use it.

16:24 Here's the way to write it more formally, and this is what's called the Kalman filter. I can basically say that my prediction x hat, so I can go back here, look at this expression here. I'm going to go ahead and add and subtract x naught here, manipulate some algebra, and here's what you get. My prediction, x bar, my improved assimilated prediction, is what my model tells me to predict, plus K, which is the Kalman term, y minus x naught. Y is my measurement, x naught is my model prediction, y is my measurement, and what I do is take the difference this times K, and K is given by this quantity here, which is less than or equal to one, and what it is is the variance of my model squared over the variance of my model squared plus the variance of my measurement squared.

17:19 One really important thing to note here, in this common filter situation, if I have a perfect measurement, in other words, the sensor I have has no variance. It's a perfect measurement. In other words, sigma y is zero. There is no variance in the measurement. If this is zero, then you have sigma naught squared over sigma naught squared. This is one. And if this is one, then you have x naught plus one y minus x naught. The x naughts cancel, and all you have is y.

17:50 In other words, what it would say, if I have a perfect measurement, my assimilated solution is simply take on the value of the measurement, because it is perfect. Okay? So, this makes a lot of sense. Hopefully, that makes a lot of sense to you, okay? Versus also, what if my measurements were really poor? I my sensors are awful. In other words, sigma y squared is really big. So, if the variance of the measurement is really big, you have these numbers divided by something big, k is going towards zero, let's say, if sigma y squared gets very large. Then, if this goes towards zero, then you take on the value instead of your model. So, in these two limits, this completely makes sense. If you have terrible, terrible measurements, use your model.

18:38 If you have really great measurements, you're going to really weight that measurement much more strongly than, let's say, your model. So, it makes a great deal of intuitive sense, and of course, what you hope to do is actually get good estimates of these variances in practice. Okay. K, this common filter, is called the innovation. In other words, I take my model, I make prediction, and I update my model with this term. That is the innovation that we would have.

19:09 So, as you can see from this, common filtering gives a pathway towards building better models by using the model we have, but then pinning our predictions down to reality, which are the reality being measurements that we get from sensors into the system. So, anyway, data simulation in a little nutshell. That's a very basic example. We'll build more here in a moment. >>