Optimal Spatial Prediction with Kriging

Robert Nishihara · April 21, 2013

Blog   Machine Learning   Statistics

Suppose we are modeling a spatial process (for instance, the amount of rainfall around the world, the distribution of natural resources, or the population density of an endangered species). We've measured the latent function $Z$ at some locations ${\mathbf{s}}_1,\dots ,{\mathbf{s}}_N$, and we'd like to predict the function's value at some new location ${\mathbf{s}}_0$. Kriging is a technique for extrapolating our measurements to arbitrary locations. For an in-depth discussion, see Cressie and Wikle (2011). Here I derive Kriging in a simplified case.

I will assume that $Z$ is an intrinsically stationary process. In other words, there exists some semivariogram $\gamma (\mathbf{h})$ such that

$$\text{var}[Z(\mathbf{s}+\mathbf{h})-Z(\mathbf{s})]=2\gamma (\mathbf{h}).$$

Furthermore, I will assume that the process is isotropic, (i.e. that $\gamma (\mathbf{h})$ is a function only of $||h||$). As Andy described here, the existence of a covariance function implies intrinsic stationarity. In addition, I will assume that the process has a constant mean, $\mathbb{E}[Z(\mathbf{s})]=\mu$. We would like to estimate $Z(\mathbf{s})$ with a linear combination of our current observations. Our estimator will be

$$\hat{Z}=\sum _{n=1}^N{\lambda }_{n}Z({\mathbf{s}}_n),$$

where the weights ${\lambda }_n$ can be positive or negative. We further require that $\sum _n{\lambda }_n=1$ so that our estimate is unbiased. We would like to choose the weights ${\lambda }_n$ so as to minimize the mean-squared predictive error

$$\text{MSPE}({\lambda }_1,\dots ,{\lambda }_N)=\mathbb{E}[(\hat{Z}-Z(\mathbf{s}){)}^2].$$

Let ${\gamma }_{nm}$ denote $\gamma ({\mathbf{s}}_n-{\mathbf{s}}_m)$. Expanding the expression for the mean-squared predictive error, we get

$$\sum _{n,m}{\lambda }_n{\lambda }_m\mathbb{E}[Z({\mathbf{s}}_n)Z({\mathbf{s}}_m)]+\mathbb{E}[Z({\mathbf{s}}_0{)}^2]-2\sum _n{\lambda }_n\mathbb{E}[Z({\mathbf{s}}_n)Z({\mathbf{s}}_0)].$$

Adding and subtracting $\sum _n{\lambda }_n\mathbb{E}[Z({\mathbf{s}}_n{)}^2]$, this expression breaks into $A+B$, where

$$\begin{array}{rcl}A& =& -\sum _n{\lambda }_n\mathbb{E}[Z({\mathbf{s}}_n{)}^2]+\sum _{n,m}{\lambda }_n{\lambda }_m\mathbb{E}[Z({\mathbf{s}}_n)Z({\mathbf{s}}_m)]\\ & =& -\frac{1}{2}\sum _{n,m}{\lambda }_n{\lambda }_m\mathbb{E}[(Z({\mathbf{s}}_n)-Z({\mathbf{s}}_m){)}^2]\\ & =& -\sum _{n,m}{\lambda }_n{\lambda }_m{\gamma }_{nm},\end{array}$$

and

$$\begin{array}{rcl}B& =& \mathbb{E}[Z({\mathbf{s}}_0{)}^2]-2\sum _n{\lambda }_n\mathbb{E}[Z({\mathbf{s}}_n)Z({\mathbf{s}}_0)]+\sum _n{\lambda }_n\mathbb{E}[Z({\mathbf{s}}_n{)}^2]\\ & =& \sum _n{\lambda }_n\mathbb{E}[(Z({\mathbf{s}}_0)-Z({\mathbf{s}}_n){)}^2]\\ & =& 2\sum _n{\lambda }_n{\gamma }_{0n}.\end{array}$$

We minimize the quantity $A+B$ subject to the constraint $\sum _n{\lambda }_n=1$ using Lagrange multipliers. To simplify the notation, define the matrix $\mathrm{\Gamma }$ by ${\mathrm{\Gamma }}_{nm}={\gamma }_{nm}$, the vector $\mathit{\lambda }=({\lambda }_1,\dots ,{\lambda }_N{)}^{\mathsf{T}}$, the vector ${\mathit{\gamma }}_0=({\gamma }_{01},\dots ,{\gamma }_{0N}{)}^{\mathsf{T}}$, and the vector $\mathbf{1}=(1,\dots ,1{)}^{\mathsf{T}}$. Then the mean-squared predictive error is given by

$$-{\mathit{\lambda }}^{\mathsf{T}}\mathrm{\Gamma }\mathit{\lambda }+2{\mathit{\lambda }}^{\mathsf{T}}{\mathit{\gamma }}_0.$$

Incorporating the Lagrange multiplier constraint, we have the quantity

$$\mathrm{\Phi }(\mathit{\lambda },\alpha )=-{\mathit{\lambda }}^{\mathsf{T}}\mathrm{\Gamma }\mathit{\lambda }+2{\mathit{\lambda }}^{\mathsf{T}}{\mathit{\gamma }}_0-\alpha ({\mathbf{1}}^{\mathsf{T}}\mathit{\lambda }-1).$$

Differentiating $\mathrm{\Phi }$ with respect to $\alpha$ gives back our constraint. Differentiating with respect to each ${\lambda }_n$ and concatenating the resulting equations into matrix form gives

$$\begin{array}{rcl}& & -2\mathrm{\Gamma }\mathit{\lambda }+2{\mathit{\gamma }}_0=\alpha \mathbf{1}\\ & \iff & \mathit{\lambda }={\mathrm{\Gamma }}^{-1}({\mathit{\gamma }}_0-\alpha \mathbf{1}/2).\end{array}$$

Incorporating the constraint gives

$$\begin{array}{rcl}1& =& {\mathbf{1}}^{\mathsf{T}}\mathit{\lambda }\\ & =& {\mathbf{1}}^{\mathsf{T}}{\mathrm{\Gamma }}^{-1}({\mathit{\gamma }}_0-\alpha \mathbf{1}/2).\end{array}$$

Solving for $\alpha$ and plugging this back into our formula for $\mathit{\lambda }$, we find that

$$\mathit{\lambda }={\mathrm{\Gamma }}^{-1}({\mathit{\gamma }}_0-\frac{{\mathbf{1}}^{\mathsf{T}}{\mathrm{\Gamma }}^{-1}{\mathit{\gamma }}_0-1}{{\mathbf{1}}^{\mathsf{T}}{\mathrm{\Gamma }}^{-1}\mathbf{1}}\mathbf{1}).$$

This gives us our optimal Kriging predictor.

Share: Twitter, Facebook