12  Details on kriging predictions


Download slides PDF

The PDF contains what was shown during the lecture. Some of these websites include additional text in light gray.

12.1 Prediction problem formulation

Observations \(\mathbf{y}^\mathrm{T} = (y_{1}, \ldots, y_{n})\) available for a set of \(n\) locations \(\mathbf{x}_{i}\)

Consider \(\mathbf{y}\) as a realization of the random variable \(\mathbf{Y}^\mathrm{T} = (Y_{1}, \ldots, Y_{n})\)

 

Model: \(Y_{i} = S(\mathbf{x}_{i}) + Z_{i}\) with:

  • \(Y_{i}\): the \(i^\mathrm{th}\) datum
  • \(S(\mathbf{x}_{i})\): “signal” (the true quantity) at location \(\mathbf{x}_{i}\)
  • \(\{S(\mathbf{x}_{i})\}\): Gaussian process, parametrized by:
    • Trend: \(\mu(\mathbf{x}_{i}) = \sum_{k} d_{k}(\mathbf{x}_{i}) \beta_{k} = \mathbf{d}(\mathbf{x}_{i})^\mathrm{T} \mathbf{\beta}\)
    • Covariance function \(\gamma(\mathbf{h}; \mathbf{\theta})\) or variogram \(V(\mathbf{h}; \mathbf{\theta})\)
  • \(Z_{i}\): independent, identically distributed (iid) Gaussian measurement error with variance \(\tau^2\)

Predictions: Let’s say \(\hat{\mathbf{S}}\) is the prediction of \(\mathbf{S}^\mathrm{T} = (S(\mathbf{x}_{1}^\prime), \ldots, S(\mathbf{x}_{m}^\prime))\) for a set of \(m\) locations \(\mathbf{x}_{j}^\prime\) without data.

\(\hat{\mathbf{S}}\) is computed from \(\mathbf{Y}\), therefore \(\hat{\mathbf{S}} = \hat{\mathbf{S}}(\mathbf{Y})\).

12.2 Kriging prediction at new point location

Ordinary kriging – with a global mean

Ordinary kriging (OK) assumes an unknown but constant global mean (just an intercept, no coefficients in the trend model). By requiring the kriging weights to sum to one, any constant component – including the unknown mean – is reproduced automatically. The prediction can therefore be written simply as a weighted sum of the observations.

Puncutal kriging: For same spatial support – prediction entities are assumed to have same extension as observations (i.e. sampling area, sensor size)

Prediction of signal \(S(\mathbf{x}_{0})\) at location \(\mathbf{x}_{0}\) without measurement

\[ \hat{S}(\mathbf{x}_{0}) = \sum_{i=1}^n \kappa_{i}(\mathbf{x}_{0})\, y(\mathbf{x}_{i}) \]

with kriging weights \(\kappa_{i}(\mathbf{x}_{0})\). To ensure unbiased estimates, weights are made to sum to 1:

\[ \sum_{i=i}^N \kappa_i(x_0) = 1 \]  

Expected difference \(\mathbb{E}[S(\mathbf{x}_{0}) - \hat{S}(\mathbf{x}_{0})] = 0\).

Estimate of kriging variances

Estimate of the variance of the prediction error:

\[ var[\hat{S}(\mathbf{x}_{0})] = \mathbb{E}\left[\left\lbrace S(\mathbf{x}_{0}) - \hat{S}(\mathbf{x}_{0}) \right\rbrace\right] \]

\[ = 2\sum_{i=1}^N \kappa_{i}(\mathbf{x}_{0}) V(\mathbf{x_i}, \mathbf{x_0}) - \sum_{i=1}^N \sum_{j=1}^N \kappa_{i}(\mathbf{x}_{0}) \kappa_{j}(\mathbf{x}_{0}) V(\mathbf{x_i}, \mathbf{x_j}) \] with \(V(\mathbf{x_i}, \mathbf{x_0})\) being the semivariance of \(S\) between sampling point \(x_i\) and the target prediction point \(\mathbf{x_0}\) and \(V(\mathbf{x_i}, \mathbf{x_j})\) the semivariance between the \(i\)th and the \(j\)th sampling points.

Kriging variance quantifies the uncertainty of the prediction at each new location. It is determined by the spatial covariance/variogram model and by the locations of the observations relative to the prediction location. The formula accounts for how much information the surrounding observations provide and how strongly they are spatially correlated.

If the observations are moved to a different configuration, the kriging variance changes – even if exactly the same measured values are kept. The uncertainty is therefore a property of the sampling geometry and the assumed spatial dependence.

Find kriging weights

Kriging weights are calculated by solving the system of equations that minimizes the prediction error variance subject to the constraints of the variogram model.

Find weights \(\kappa_{i}(\mathbf{x}_{0})\) that minimize kriging variances and sum to 1 by solving \(N + 1\) equations:

\[ \sum_{i=1}^N \kappa_{i}(\mathbf{x}_{0}) V(\mathbf{x_i}, \mathbf{x_j}) + \psi(\mathbf{x_0}) = V(\mathbf{x_j}, \mathbf{x_0}) \quad \text{for all} ~~j \]

\[ \sum_{i=i}^N \kappa_i(x_0) = 1 \]

The quantity \(\psi(\mathbf{x_0})\) is a Lagrange multiplier introduced to achieve minimization. Computed as the inverse of the semivariance matrix multiplied by the vector of semivariances to the target point.

The kriging weights depend on the prediction location. For each new location, the kriging equations are solved using the distances (or covariances) between that location and all observations. This gives a new set of weights for every prediction location, which are then used to calculate the kriging prediction and its uncertainty.

12.3 Universal/external drift kriging predictions

Universal kriging: often referred if trend is modelled by coordinates

External-drift kriging: trend is modelled by spatial covariates


Evaluating \(\hat{\mathbf{S}}_{\text{opt}}\) requires a fully specified weakly stationary model:

  1. Structure of trend function is known

  2. Regression coefficients \(\mathbf{\beta}\) are known

  3. Type of parametric covariance (variogram) function is known

  4. Parameters \(\mathbf{\theta}\) and \(\tau^2\) of the covariance function are known

Relaxed assumptions: Only 1, 3, and 4 are assumed to be known, while \(\mathbf{\beta}\) is implicitly estimated from the data using generalized least squares (GLS).

UK/EDK predictions

Therefore, we use the universal kriging (UK) or external drift kriging (EDK) plug-in predictor: \[ \hat{\mathbf{S}}_{\text{k}} = \mathbf{X}_{\mathbf{S}} \hat{\mathbf{\beta}}_{\text{GLS}} + \mathbf{\Lambda} \left( \mathbf{y} - \mathbf{X}_{\mathbf{Y}} \hat{\mathbf{\beta}}_{\text{GLS}} \right) \]
with: \[ \mathbf{\Lambda} = \mathbf{\Sigma}_{\mathbf{S}\mathbf{Y}} \mathbf{\Gamma}_{\mathbf{Y}\mathbf{Y}}^{-1} \]


Computing universal kriging predictor requires:

  1. known structure of trend function

  2. known structure and parameters of variogram

“Plug-in” predictor: uncertainty of variogram is ignored when computing predictions

Universal kriging (UK) can be viewed as trend prediction plus spatial prediction of the residuals. First, the trend is estimated using GLS. The residuals are then predicted at the new location using their spatial covariance with the observations. The covariance weights are \(\Lambda=\Sigma_{SY}\Gamma_{YY}^{-1}\). Finally, the predicted trend and predicted residual are added.

UK follows the idea of ordinary kriging, but assumes a known trend model. So it is closer to simple kriging (with a known mean, not discussed in this course).

Unlike ordinary kriging, the unbiasedness constraints must reproduce the specified trend functions rather than only a constant mean (constraint to sum to 1). In the plug-in formulation, the trend is estimated first by GLS and the remaining spatial residual is then predicted from its covariance with the observations through \(\Lambda\).

Matrix

\[ \Gamma_{YY}=\operatorname{Cov}(Y,Y^\top) \]

is the covariance matrix among the observations. Its \(ij\)-th element is the covariance between observations at locations \(s_i\) and \(s_j\).

And

\[ \Sigma_{SY}=\operatorname{Cov}(S_k,Y^\top) \]

is the cross-covariance between the quantity being predicted at \(s_k\) and the observations.

Therefore,

\[ \Lambda=\Sigma_{SY}\Gamma_{YY}^{-1} \]

gives the covariance-based weights used to predict the spatial residual at \(s_k\).