Spatial Values-Spatial Prediction

Hi, I'm Serna, a consultant at Wolfram Research in the Probability and Statistics group. Like OSI said, there's two of us in the group and we've been here for many years and this is my talk on spatial values and spatial prediction. This is a similar talk to one I gave last year. I've tried to emphasize some quality of life improvements that I've added since last release and also showcase some nice examples that are upcoming. Functionality in different areas of the of the language. So to begin, spatial data sets consisting of a set of measured values at specific locations are becoming increasingly important. Some examples include temperature, elevation, concentration of minerals, etc. We're going to have a look at existing Wolfram language functionality to perform estimation of missing values in a region of interest. Modeling an underlying quantity Traditionally this procedure has been known as kriging. So first example, let's have a look at a spatial data set that is a set of continuous quantity with large gaps. Ocean readings over the contiguous United States. This is a perfect example for spatial estimates because it's an estimator that can combine both in large scale trends and local variation. That means that it can fill in the gaps where we're very far from points in the data set, while at the same time capture the structure when we're closer to a point in the data set. So let's just have a look at the first. Have a look at the points in the data set and we see that there's a large gap. This is presumably because the data collecting this, these satellites collecting this data, the orbits don't close on themselves, so there's some gaps in their measurements. Let's create a spatial estimator function with a large scale quadratic trend. And this trend is specified by the option spatial trend function. And we use the integer 2 to specify the order of the fitted polynomial that we're going to use for the trend. The estimator function captures. The summary box here captures many things like the spatial noise level, the order of the polynomial, the type of the variogram, these I'm going to talk about these in a little bit and the representation of the symbolic trend. We can evaluate the spatial predictor function at uniform random locations in the observation region and use these locations, these evaluated locations, to compare the plots between the data set on the left, the original data set on the left, and the plot of the estimator on the right. And we can see that indeed the estimating is doing a remarkable job of literally filling in the gaps. We're capturing a lot of Structure that was coming both trends and more local variations that were coming from the points in the data set. We can also for instance, use this estimator to just compute an estimated value in a particular location. For instance, London, Ed. Not London, Boston. The position of Boston. And we get the reading of the ozone in whatever units this was in to begin with. Okay. As a warm up I'm going to use the familiar model of ordinary least squares. And the model is usually as follows. We have a function, a parametric model F, some locations x, some values Y, both in x X, a part of the data set. And we're trying to estimate the parameters of F. And we assume there's some unknown noise delta that we don't control, where here delta is identically independently distributed zero mean Gaussian random variable, and F is the function that we're trying to estimate. Once F is estimated, we have the estimated, the parameters are estimated, we call it F hat, and we can infer the value of the model Y hat at a particular location. But just computing the value of f hat at that location. Let's have a look at an example. For instance, I've chosen this data so that it is plausible and reasonable to assume that this is a some kind of polynomial with some noise introduced to it. And we can use the function linear model fit, which has been in the system for a long time, to fit a cubic polynomial to it. Here we're giving the basis functions and it will estimate form linear squares, least squares on the system and gets a fitted model. We can plot this model using the function accessories and we see that indeed it's doing what it's supposed to be doing. It captures a reasonable polynomial very well. A cubic seems to be a reasonable model for this data. Let's have a look now at the model part of the of the estimator. And in particular let's have a look at the DF in particular here. This model in spatial estimate is essentially the same as in least squares. In fact, we're using generalized least squares to minimize. But instead of using F hat as our estimator, we're using a linear combination of the Y values with some weights lambda that we are going to estimate and then depend on the point of interest and are chosen to minimize the expected squared difference between the estimator and the model itself. This is called an unbiased minimum variance linear estimator. Unbiased because it minimizes the it in minute. Okay, one second. It's unbiased because it has the same mean as the as the original model and it minimizes the variance because it's minimizing this quantity. And it's a linear estimator, because this is a linear estimator. Let's have a look at the previous example and use spatial estimates on the data and the values that we had before and plot it alongside with the plot that we had before. And we can see that our estimator in this case is an interpolator and it's going through every single data point as opposed to the cubic which is just going close to a data points. Let's now have a look at the noise part of the model component and we can further split the 0 mean random component delta into a correlated component eta and a white noise identically independent distributed component epsilon which has zero covariance as opposed to the e term component which will have covariance non vanishing covariance between points that are in different locations. The variance of epsilon is specified by the option spatial noise level and it is a measure of how noisy the data is or in other words how much our estimator should smooth things out. If the data is noisier, we are less sure of our estimator being close to the data. Let's compute different estimators with an increasing value of the noise, an increasing level of noise and up to a reasonably small value, up to an increasingly small value. And let's have. Let's just plot them and see how they look. And we can see that in fact the zero noise level is an exact interpolator. And then as the, as the. As the noise rises, we're getting closer and closer to what used to be the trend. And in fact. This what's going on the. When the trend F, which is estimated according to generalized least squares, as the level of noise goes to infinity, the spatial estimator approaches the approaches the F hat our previous generalized least squares estimator. This is. So let's compute two estimators with a very large level of the noise 10 and then compute both one with a cubic trend and a quadratic trend. And if we plot them we see that indeed we have that we can see in the red line here that the cubic trend is exactly the same as the estimated with the cubic trend is exactly the same as the cubic fits model that we had before and that indeed the estimators are just essentially the trend. We can also query the estimator for the trend using the symbolic trend accessor. This is also in the summary box and it will just give us a symbolic representation of the trend and we can use the GLS trend, which stands for generalized least squares. We can use this accessor to obtain a functional representation of the trend that we can use for plotting purposes or just further analysis. Okay, so just let's have a look at the last part of our model. The eta and which is, which is specified by the what we call the variogram gamma. And this is for a stationary random process. This gamma is a function of only the distance. The gamma is given by the expected squared difference between the random component eta at different locations. And let's look at an example of measurements of rainfall in Switzerland to see how the variogram affects the estimator. We can get this data using the resource data function and it returns a spatial point data object, which my colleague is going to talk about more in her talk in a minute. And she and the. We can use the observation region accessor to obtain the region of observation and we can feed it directly to point value plot to get a plot of the intensity and the location, the intensity of the different the amount of rainfall at the different weather stations where these measurements were taken. We can also use the annotations of the spatial point data object directly to feed it to spatial estimate to get the appropriate value to use. So we can query the annotations like this. And we see that we only have one annotation and we can use this string instead of a list of values. And this will do what we expected spatial estimate by default will automatically estimate the variogram. Here we see that it's estimated a qubit variogram with this level of noise and these parameters that are meaningless until you look at a plot. And when we do so we see that the plot and for this case, I mean I did choose the example looks very nice. The level of the variogram is very close to the empirical variagram in this case, and it's very reasonable. We can obtain this empirical variagram using the function binned variogram list with the same input as before. And when we plot it and this returns an awaited data object. And when we plot it we just get what we have here in the visualization of variagram the points. And by default a suitable model chosen, the cubic model has some very nice properties, but we can also manually specify a varagar model to fit. And this is done by using an estimated varicam model. And in this case I'm going to choose the linear model. The reason I'm choosing this one is because it looks okay, but worse than what we had before, and I'm doing this on purpose. We can see that it's reasonable, it could be linear, but it's clearly visually worse than what we, than the model we had before. And also quantitatively we can compute the, we can access, query the weighted sum of residuals and see that indeed in the case for the linear one, it's actually much worse. So now let's just plot these two different estimators. Let's compute the estimator for the linear case, for the linear variogram and plot them side by side. And we can see that the cubic one in fact captures more detail and is, you know, you can see a lot more structure in the data. So choosing the right variable matters, but we also try and do it for you as best as possible. These models must have the property that that amounts to the fact that the covariance of the points has to be positive, definite. And this is a hard problem to come up with these functions, a hard open research problem. So there's some popular models that we just provided for you so that you don't have to do it, code them by hand. Okay, let's revisit spatial noise level again. I wanted to point out some of the things that spatial estimate is doing under the hood and some of the subtleties with different variables and spatial noise level. Let's have a look at the Gaussian variagram, which is a popular variagram for reasons, but it also has some difficulties in that the flat slope at the origin causes it to suffer from conditioning problems if the data points are too close together. This means that the linear system that the spatial estimate has to solve will be ill conditioned. And if you don't regularize the system appropriately or raise precision to ensure that you do not lose precision, you're going to get gibberish. So let's compute two estimators, both. Well, let's first specify two variograms, one with no, with zero noise and one with where we don't specify the noise. And this means that the system is going to do it for us and the system is going to try and make sure that the system is not ill conditioned and that we're not losing precision needlessly. So we compute these two estimators with these two different variagrams. And we're also going to set precision to 16 digits and we're going to do the same with this time at arbitrary precision. This means that what the system's going to do is it's going to raise precision under the hood when it needs to. To make sure that we don't lose any precision. It's going to compute internally, compute the conditioning number of the system. And let's just have a look at the different plots. And we can see that here on the top right we have the regularized automatic noise level. And the plot of the estimator is very smooth and very reasonable. It goes close to the data, but it's nice. On the left we have the one where we haven't regularized. And we also have fixed precision. And it's gibberish. We see it's the values are very, very far away from the data, which are of order one. It's nonsense. On the bottom right we have the arbitrary precision case, but we've allowed the system to estimate the noise level. And it's chosen a noise level not by regularizing, but by looking at the plausible noise of the data. And it's decided that by fitting the variogram, essentially when it fits the variogram, it estimates the noise as well. And the noise, the noise that it's come up with is a reasonable regularizer of the problem as well. So we also get a nice plot and for the case where we choose where we're allowed to raise precision, but we keep the noise at zero, we can see. I'm going to close up here below this is the same plot and this is a close up. We see that we're getting an estimator that does interpolate, but it's kind of useless. It oscillates wildly outside of the far away from the data sets. It goes up to infinity very fast. And this is not a very useful thing. Okay, how am I doing in time? Very close. Okay, almost done. Okay, so let's look at a nice example now. Oops. This initialization was supposed to not be here. Okay, I'm sorry, I can't minimize it now. Yes I can. Okay, the. Let's look at. Let's use this upcoming functionality of astroposition that my Jose I said they talked about or is going to talk about in the in and a talk. And this we're going to use entity class to obtain a list of the 10 brightest stars in the night sky. And these stars are. We can use astral position to obtain their coordinates in the ICRS inertial frame. And we take the first two returns a list of triplets. The first two are the angles in spherical coordinates and the third one is the distance from the coordinate origin. The angles are in essentially latitude and longitude, although you can. They're actually called Right ascension and declination. We use apparent magnitude to capture the brightness from Earth. So it's not the intrinsic brightness, but how with apparent brightness. Apparent brightness is a. Apparent magnitude is a function of both how bright intrinsically bright it is, and the distance to us. And the brightness is proportional to the negative to the exponential of the negative of the apparent magnitude. So we use the apparent magnitude accessor to query the magnitude. And if we plot it, if we plot the brightness, we see that indeed it's a decreasing quantity, reasonably decreasing quantity. And we're going to compute a spatial estimator. We need to be able to compute the distances in the unit sphere. So we're going to use the geodesic here in the unit sphere in terms of the angles of the latitude and longitude. And this is the appropriate expression. You can use the Haversine formula, that is better condition. But this works fine for this example. So we compute the spatial estimator using the brightness, the spherical geodesic function, and we use a quadratic trend. And then we compute some locations in the celestial sphere to interpolate, to compute our estimator on. And we can then use a point, produce a plot of these, both of the original stars and their brightness and the estimator that we just computed. And this allows us, we can, we can wrap this texture on the sphere and we can have an estimated like here, this star is presumably Sirius is much brighter than the others. And we can produce an interpolator, an estimator of the brightness of the night sky. There's a here, there's, you know, the texture doesn't quite wrap all the way around in the singularities in the coordinate system. That is not a big deal unless you want it to be pretty. To conclude, this functionality has been in the system since 13.0 and we hope it has been useful to you guys. We have been making some quality of life improvements. We hope you keep using it and you let us know what features you're interested in. Thank you very much.