8
$\begingroup$

I have followed tutorials online for spatial kriging with both geoR and gstat (and also automap). I can perform spatial kriging and I understand the main concepts behind it. I know how to build a semivariogram, how to fit a model to it and how to perform ordinary kriging.

What I don't understand is how the weights of the surrounding measured values are determined. I know they derive from the semivariogram and depend on the distance from the prediction location and on the spatial arrangement of the measured points. But how?

Could anyone please make an ordinary kriging (non-bayesian) model with 3 random measured points and 1 prediction location? It would be enlightening.

$\endgroup$
  • 1
    $\begingroup$ just for curiosity, why don't you want to see the Bayesian answer? It makes things much simpler when you deal with Gaussian Processes. $\endgroup$ – DeltaIV Feb 7 '18 at 10:13
  • $\begingroup$ @DeltaIV because first I want to learn the frequentist way. Bayesian statistics are still cloudy for me $\endgroup$ – Pigna Feb 7 '18 at 10:47
  • 1
    $\begingroup$ "What I don't understand is how the weights of the surrounding measured values are determined.". In case anyone is interested, I posted an answer in GIS SE with an example about how to calculate them (gis.stackexchange.com/questions/270274/…). But the answer here is great already! $\endgroup$ – Andre Silva Nov 5 '18 at 14:39
11
+50
$\begingroup$

Apart from this answer, there are also some nice additional answers to a similar question on gis.stackexchange.com

First I'll describe ordinary kriging with three points mathematically. Assume we have an intrinsically stationary random field.

Ordinary Kriging

We're trying to predict the value $Z(x_0)$ using the known values $Z=(Z(x_1),Z(x_2),Z(x_3))$ The prediction we want is of the form $$\hat Z(x_0) = \lambda^T Z$$ where $\lambda = (\lambda_1,\lambda_2,\lambda_3)$ are the interpolation weights. We assume a constant mean value $\mu$. In order to obtain an unbiased result, we fix $\lambda_1 + \lambda_2 + \lambda_3 = 1$. We then obtain the following problem: $$\text{min} \; E(Z(X_0) - \lambda^T Z)^2 \quad \text{s.t.}\;\; \lambda^T \mathbf{1} = 1.$$ Using the Lagrange multiplier method, we obtain the equations: $$\sum^3_{j=1} \lambda_j \gamma(x_i - x_j) + m = \gamma(x_i - x_0),\;\; i=1,2,3,$$ $$\sum^3_{j=1} \lambda_j =1 ,$$ where $m$ is the lagrange multiplier and $\gamma$ is the (semi)variogram. From this, we can observe a couple of things:

  • The weights do not depend on the mean value $\mu$.
  • The weights do not depend on the values of $Z$ at all. Only on the coordinates (in the isotropic case on the distance only)
  • Each weight depends on location of all the other points.

The precise behaviour of the weights is difficult to see just from the equation, but one can very roughly say:

  • The further the point is from $x_0$, the lower its weight is ("further" with respect to other points).
  • However, being close to other points also lowers the weight.
  • The result is very dependent on the shape, range, and, in particular, the nugget effect of the variogram. It would be quite illuminating to consider kriging on $\mathbb R$ with only two points and see how the result changes with different variogram settings.

I will however focus on the location of points in a plane. I wrote this little R function that takes in points from $[0,1]^2$ and plots the kriging weights (for exponential covariance function with zero nugget).

library(geoR)

# Plots prediction weights for kriging in the window [0,1]x[0,1] with the prediction point (0.5,0.5)
drawWeights <- function(x,y){
 df <- data.frame(x=x,y=y, values = rep(1,length(x)))
  data <- as.geodata(df, coords.col = 1:2, data.col = 3)
  
  wls <- variofit(bin1,ini=c(1,0.5),fix.nugget=T)
  weights <- round(as.numeric(krweights(data$coords,c(0.5,0.5),krige.control(obj.mod=wls, type="ok"))),3)
  
  plot(data$coords, xlim=c(0,1),  ylim=c(0,1))
  segments(rep(0.5,length(x)), rep(0.5,length(x)),x, y, lty=3 )
  text((x+0.5)/2,(y+0.5)/2,labels=weights)
}

You can play with it using spatstat's clickppp function:

library(spatstat)
points <- clickppp()
drawWeights(points$x,points$y)

Here are a couple of examples

Points equidistant from $x_0$ and from each other

deg <- seq(0,2*pi,length.out=4)
deg <- head(deg,length(deg)-1)
x <- 0.5*as.numeric(lapply(deg, cos)) + 0.5
y <- 0.5*as.numeric(lapply(deg, sin)) + 0.5
drawWeights(x,y)

enter image description here

Points close to each other will share the weights

deg <- c(0,0.1,pi)
x <- 0.5*as.numeric(lapply(deg, cos)) + 0.5
y <- 0.5*as.numeric(lapply(deg, sin)) + 0.5
drawWeights(x,y)

enter image description here

Nearby point "stealing" the weights

deg <- seq(0,2*pi,length.out=4)
deg <- head(deg,length(deg)-1)
x <- c(0.6,0.5*as.numeric(lapply(deg, cos)) + 0.5)
y <- c(0.6,0.5*as.numeric(lapply(deg, sin)) + 0.5)
drawWeights(x,y)

https://i.imgur.com/MeuPvFT.png

It is possible to get negative weights

enter image description here

Hope this gives you a feel for how the weights work.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.