We can think of a Gaussian process as a collection of random variables that follow a multivariate Gaussian distribution, as in \[y(s) = \beta_0 + e(s), e(s) \sim \mathcal{N}(0,\Sigma)\] with \(\Sigma_{ij}=Cov(e(s_i),e(s_j))\) (which is a variance if \(i=j\)). The index \(s\) indicates that \(y(s)\) might be a spatial, temporal, or spatiotemporal process. \(\beta_0\) can be estimated by generalized least squares, but the interest usually lies in predicting (estimating) \(y(s_0)\) at an unobserved location \(s_0\) from \(n\) observed values \(y(s_i), i = 1,...,n\). Given known \(\beta_0\) and \(\Sigma\), and with \(\Sigma_0\) the covariance between scalar \(y(s_0)\) and vector \(y(s)\), we get \[\hat{y}(s_0)=\beta_0 + \Sigma_0 \Sigma^{-1}(y(s)-\beta_0)\] which reduces to \(\beta_0\) when \(y(s_0)\) is uncorrelated with the other \(y(s)\) (and \(\Sigma_0\) is a vector with zeroes). This predictor is the best linear predictor, Gaussian process regression, and in spatial context also known as simple kriging.
In more realistic situations, \(\beta_0\) is not known (ordinary kriging) or a linear model \(X(s)\beta\) with unknown parameter vector \(\beta\) is assumed for the mean function (universal kriging). In that case, the generalized least squares (GLS) estimate for \(\beta=(X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1}y(s)\) is substituted, giving \[\hat{y}(s_0)=x_0\hat{\beta}+\Sigma_0\Sigma^{-1}(y(s)-X\hat{\beta})\] with \(x_0\) the vector with predictor values corresponding to \(s_0\). In practice, when e.g. spatially predicting values over a regular grid covering the target region, for each \(s_0\) only \(x_0\) and \(\Sigma_0\) change, all other parts remain constant.
The theory of GP is pretty straightforward, the practice a bit harder:
\[ \newcommand{\E}{{\rm E}} % E expectation operator \newcommand{\Var}{{\rm Var}} % Var variance operator \newcommand{\Cov}{{\rm Cov}} % Cov covariance operator \newcommand{\Cor}{{\rm Corr}} \] ### Unknown, constant mean Suppose the mean is constant, but not known. This is the most simple realistic scenario. We can estimate it from the data, taking into account their covariance (i.e., using weighted averaging, or GLS):
\[\hat{\beta}_0 = ({\bf 1}'\Sigma^{-1}{\bf 1})^{-1} {\bf 1}'\Sigma^{-1}y(s)\] with \({\bf 1}\) a conforming vector with ones, and substitute this mean in the SK prediction equations: BLUP/Ordinary kriging: \[\hat{y}(s_0) = \hat{\beta}_0 + v'V^{-1} (y(s)-\hat{\beta}_0)\] Which has prediction error variance \[\sigma^2(s_0) = \sigma^2_0 - \Sigma_0'\Sigma^{-1}\Sigma_0 + Q\] with \(\sigma^2_0\) = var(y(s))$ and \(Q = (1 - {\bf 1}'\Sigma^{-1}\Sigma_0)'({\bf 1}'\Sigma^{-1}{\bf 1})^{-1}(1 - {\bf 1}'\Sigma^{-1}\Sigma_0)\)
Given prediction location \(s_0\), and data locations \(s_1\) and \(s_2\), we need: \(\Var(y(s_0))\), \(\Var(y(s_1))\), \(\Var(y(s_2))\), \(\Cov(y(s_0),y(s_1))\), \(\Cov(y(s_0),y(s_2))\), \(\Cov(y(s_1),y(s_2))\).
How to get these covariances?
The solution is to assume stationarity.
Stationarity of the
Second order stationarity: \(\Cov(y(s),y(s+h)) = C(h)\)
which implies: \(\Cov(y(s),y(s)) = \Var(y(s))= C(0)\)
The function \(C(h)\) is called the covariogram of the random function \(y(s)\)
Covariance: \(\Cov(y(s),y(s+h)) = C(h) = \E[(y(s)-m)(y(s+h)-m)]\)
Semivariance: \(\gamma(h) = \frac{1}{2} \E[(y(s)-y(s+h))^2]\)
\(\E[(y(s)-y(s+h))^2] = \E[(y(s))^2 + (y(s+h))^2 -2y(s)y(s+h)]\)
\(\E[(y(s)-y(s+h))^2] = \E[(y(s))^2] + \E[(y(s+h))^2] - 2\E[y(s)y(s+h)] = 2\Var(y(s)) - 2\Cov(y(s),y(s+h)) = 2C(0)-2C(h)\)
\(\gamma(h) = C(0)-C(h)\)
\(\gamma(h)\) is the semivariogram of \(y(s)\).
library(sf)Linking to GEOS 3.12.2, GDAL 3.11.4, PROJ 9.4.1; sf_use_s2() is TRUEdata(meuse, package = "sp")
meuse = st_as_sf(meuse, coords = c("x", "y"))
library(gstat)
v = variogram(log(zinc) ~ 1, meuse)
v.fit = fit.variogram(v, vgm(1, "Sph", 900, 1))
plot(v, v.fit)
average squared differences: \[\hat{\gamma}(\tilde{h})=\frac{1}{2N_h}\sum_{i=1}^{N_h}(y(s_i)-y(s_i+h))^2 \ \ h \in \tilde{h}\]
plot(v, v.fit)
v.fit model psill range
1 Nug 0.05066243 0.0000
2 Sph 0.59060780 897.0209vgm(psill = 0.6, model = "Sph", range = 900, nugget = 0.06) model psill range
1 Nug 0.06 0
2 Sph 0.60 900or simpler, with un-named arguments:
vgm(0.6, "Sph", 900, 0.06) model psill range
1 Nug 0.06 0
2 Sph 0.60 900Covariance: \(\Cov(y(s),y(s+h)) = C(h) = \E[(y(s)-\beta_0)(y(s+h)-\beta_0)]\)
Semivariance: \(\gamma(h) = \frac{1}{2} \E[(y(s)-y(s+h))^2]\)
\[\gamma(h)=C(0)-C(h)\]
$$ % E expectation operator % Var variance operator % Cov covariance operator
For this, we need to know how the mean varies. Suppose we model this as a linear regression model in \(p\) known predictors: \[y(s_i) = \sum_{j=0}^p \beta_j X_j(s_i) + e(s_i)\] \[y(s) = \sum_{j=0}^p \beta_j X_j(s) + e(s) = X(s)\beta + e(s)\] with \(X(s)\) the matrix with predictors, and row \(i\) and column \(j\) containing \(X_j(s_i)\), and with \(\beta = (\beta_0,...\beta_p)\). Usually, the first column of \(X\) contains zeroes in which case \(\beta_0\) is an intercept.
Predictor: \[\hat{y}(s_0) = x(s_0)\hat{\beta} + v'V^{-1} (y-X\hat{\beta}) \] with \(x(s_0) = (X_0(s_0),...,X_p(s_0))\) and \(\hat{\beta} = (X'V^{-1}X)^{-1} X'V^{-1}Z\) it has prediction error variance \[\sigma^2(s_0) = \sigma^2_0 - \Sigma_0'\Sigma^{-1}\Sigma_0 + Q\] with \(Q = (x(s_0) - X'\Sigma^{-1}Sigma_0)'(X'\Sigma^{-1}X)^{-1}(x(s_0) - X'\Sigma^{-1}\Sigma_0)\)
This form has several names: external drift kriging, universal kriging or regression kriging.
Example in meuse data set: log(zinc) depending on sqrt(meuse)
Plotting them:
As opposed to the ordinary kriging model, the universal kriging model needs knowledge of the mean vector in order to estimate the semivariance (or covariance) from the residual vector: \[\hat{e}(s) = y(s) - X\hat\beta\] but how to get \(\hat\beta\) without knowing \(V\)? This is a chicken-egg problem. The simplest, but not best, solution is to plug \(\hat{\beta}_{OLS}\) in, and from the \(e_{OLS}(s)\), estimate \(V\) (i.e., the variogram of \(y(s)\))
β¦ involves errors, uncertainties
If \(y\) has no spatial correlation, all covariances are zero and \(\Sigma_0=0\) and \(\Sigma=\mbox{diag}(\sigma^2)\). This implies that \[\hat{Z}(s_0) = x(s_0)\hat{\beta} + \Sigma_0'\Sigma^{-1} (y(s)-X\hat{\beta}) \] with \(\hat{\beta} = (X'\Sigma^{-1}X)^{-1} X'\Sigma^{-1}y(s)\) reduces to
\[\hat{y}(s_0) = x(s_0)\hat{\beta}\] with \(\hat{\beta} = (X'X)^{-1} X'y(s)\), i.e., ordinary least squares regression prediction.
Note that
We can write \[y(s)=X\beta+e(s),~ {\rm Cov}(e(s))=\Sigma\] If we write \(e(s) = u(s)+e'\), with \({\rm Cov}(u(s))=\tilde{\Sigma}\) and IID \(e' \mathcal(N)(0, \sigma_n I)\), then \(\Sigma=\tilde{\Sigma}+\sigma_n I\) and we can write our GP as a mixed model:
\[y(s) = X\beta + u(s) + e'\] with \(u(s)\) the spatially correlated random effect (a GP), and \(e'\) a white noise, IID (independently, identically distributed) residual.
For an application of universal kriging in mapping of air quality variables, and a comparison to linear regression, see e.g. Beelen et al. (2009).