Prediction equations

This appendix assumes some familiarity with the matrix notation introduced in Section 2.7. From the universal kriging equations, ordinary kriging can be derived as the special case where $p=1$ and $F$ and $f(s_0)$ contain only ones. Ordinary least squares prediction is a special case of uncorrelated weighted least squares prediction and estimation (constant weights).

Using the matrix notation of Section 2.7, the observations $z(s)$ are represented by the model

$$Z(s) = F \beta + e(s),\hspace{.5cm}{\rm E}(e(s))=0,\hspace{.5cm}{\rm Cov}(e(s)) = V$$ (A1)

with $Z(s)=(Z(s_1 ),...,Z(s_n))'$, $F = (f_1 (s),...,f_p (s))$ with
$f_i (s) = (f_i (s_1 ) , ... , f_i (s_n))'$, and $\beta = (\beta_1 ,..., \beta_p)'$. Given this model (i.e., $V$ is known), the best linear unbiased prediction (kriging predictor) of $Z(s_0)$ is

$$\hat Z (s_0 ) = f (s_0 ) \hat \beta + v_0 ' V^{-1} (z(s) - F \hat \beta),$$ (A2)

with $f(s_0) = (f_1(s_0),...,f_p(s_0 ))$ and $v_0 = ({\rm Cov}(e(s_1 ), e(s_0 )), ... , {\rm Cov}(e(s_n ), e(s_0)))'$, and where $\hat \beta$ is the best linear unbiased estimate of $\beta$:

$$\hat \beta = (F'V^{-1} F)^{-1} F'V^{-1} z(s)$$ (A3)

The kriging prediction has prediction variance (kriging variance)

$${\rm Var}(Z(s_0 )- \hat Z (s_0 )) = \sigma^2_{Z(s_0 )} - v_0 ' V^{-1} v_0 +$$ (A4)

$$(f (s_0 ) - v_0 ' V^{-1} F) (F'V^{-1} F)^{-1} (f (s_0 ) - v_0 ' V^{-1} F)'$$

with $\sigma^2_Z (s_0 ) = {\rm Var}(Z(s_0))$. When the trend contains an intercept (which will be true in most cases), it is enough to know generalised covariances, defined as $c - \gamma(h)$, for arbitrary $c$ and with $\gamma(h)$ the variogram of $Z(s)$. Gstat chooses $c$ as the sill of the variogram. When block predictions $\hat Z (B_0 )$ are required, then for user-defined base functions the values for $f(B_0 )$ should be given as input to gstat (in the mask map(s) or in the X-columns of the data() file), whereas point-to-block and block-to-block covariances (i.e., ${\rm Cov}(e(s_i ), e(B_0 ))$ and $\sigma^2_Z (B_0 )$) are derived from the point-to-point (generalised) covariances, using Gaussian quadrature [1] or, when either nblockdiscr or area is specified, using simple integration (regular or user-specified discretization, equally weighted).