4.3 Kriging con media desconocida: kriging universal y kriging residual

Como ya se comentó, el kriging universal se basa en el siguiente modelo: \[Z(\mathbf{s}) = \sum\limits_{j=0}^{p}X_{j}(\mathbf{s})\beta_{j} + \varepsilon(\mathbf{s}),\] donde \(\boldsymbol{\beta}=(\beta_{0}, \ldots, \beta_{p})^\top \in \mathbb{R}^{p+1}\) es un vector desconocido, \(\left\{ X_{j}(\cdot):j=0, \ldots,p\right\}\) son funciones conocidas y \(\varepsilon(\cdot)\) un proceso espacial de media cero con variograma conocido \(2\gamma(\mathbf{s}_{1},\mathbf{s}_{2}) = Var(\varepsilon(\mathbf{s}_{1})-\varepsilon(\mathbf{s}_{2}))\) (aunque en la práctica se suele suponer estacionario). Supondremos también que \(X_{0} (\cdot)\equiv 1\), de esta forma además en el caso particular de \(p=0\), se corresponderá con el modelo del kriging ordinario (ver Sección 4.1) muy utilizado en la práctica. Utilizando una notación matricial podemos escribir: \[\mathbf{Z}=\mathbf{X}\boldsymbol{\beta}+\boldsymbol{\varepsilon},\] siendo \(\boldsymbol{\varepsilon}=\left( \varepsilon(\mathbf{s}_{1}), \ldots, \varepsilon(\mathbf{s}_{n} )\right)^\top\) y \(\mathbf{X}\) una matriz \(n\times (p+1)\) con \(\mathbf{X}_{ij} =X_{j-1} (\mathbf{s}_{i})\), y: \[Z(\mathbf{s}_{0})=\mathbf{x}_0^\top\boldsymbol{\beta}+\varepsilon(\mathbf{s}_{0}),\] con \(\mathbf{x}_0=\left(X_{0}(\mathbf{s}_{0}), \ldots, X_{p}(\mathbf{s}_{0})\right)^\top\).

En este caso, como un predictor lineal verifica que: \[E\left( \sum\limits_{i=1}^{n}\lambda_{i} Z(\mathbf{s}_{i}) +\lambda_{0} \right) = \boldsymbol{\lambda}^\top \mathbf{X}\boldsymbol{\beta}+\lambda_{0},\]
siendo \(\boldsymbol{\lambda}=\left( \lambda_{1}, \ldots,\lambda_{n} \right)^\top\), una condición necesaria y suficiente para que el predictor sea uniformemente insesgado, i.e. \(E(p(\mathbf{Z},\mathbf{s}_{0}))=E(Z(\mathbf{s}_{0}))=\mathbf{x}_0^\top\boldsymbol{\beta}\), \(\forall \boldsymbol{\beta}\in \mathbb{R}^{p+1}\), es que \(\lambda_{0} =0\) y: \[\begin{equation} \boldsymbol{\lambda}^\top \mathbf{X} = \mathbf{x}_0^\top. \tag{4.1} \end{equation}\] Además como \(X_{0} (\cdot)\equiv 1\), una de estas restricciones es: \[\begin{equation} \sum\limits_{i=1}^{n}\lambda_{i} = 1, \tag{4.2} \end{equation}\] que es la única condición que deben verificar los pesos en el caso del kriging ordinario.

Por tanto el predictor del kriging universal será de la forma: \[p(\mathbf{Z},\mathbf{s}_{0})=\sum\limits_{i=1}^{n}\lambda_{i} Z(\mathbf{s}_{i}),\] verificando (4.1) y tal que minimiza el MSPE. Entonces se trata de minimizar: \[\begin{equation} E\left( Z(\mathbf{s}_{0})-\sum\limits_{i=1}^{n}\lambda_{i} Z(\mathbf{s}_{i}) \right)^{2} - 2\sum\limits_{j=0}^{p}m_{j} \left( \sum\limits_{i=1}^{n} \lambda_{i} X_{j} (\mathbf{s}_{i})- X_{j}(\mathbf{s}_{0}) \right) \tag{4.3} \end{equation}\] respecto a \(\left\{ \lambda_{i} :i=1, \ldots,n\right\}\) y \(\left\{ m_{j} :j=0, \ldots,p\right\}\), multiplicadores de Lagrange que garantizan (4.1). Teniendo en cuenta que el predictor es insesgado y que los pesos verifican (4.2), entonces: \[\begin{aligned} \left( \sum\limits_{i=1}^{n}\lambda_{i} Z(\mathbf{s}_{i}) - Z(\mathbf{s}_{0})\right)^2 & = \left( \sum\limits_{i=1}^{n}\lambda_{i} \varepsilon(\mathbf{s}_{i}) - \varepsilon(\mathbf{s}_{0})\right)^2 \\ & = -\dfrac{1}{2} \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\lambda_{i}\lambda_{j} \left( \varepsilon(\mathbf{s}_{i}) - \varepsilon(\mathbf{s}_{j} )\right)^{2} + \sum\limits_{i=1}^{n}\lambda_{i} \left( \varepsilon(\mathbf{s}_{i}) - \varepsilon(\mathbf{s}_{0}) \right)^{2}, \end{aligned}\] y podemos escribir (4.3) como: \[-\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\lambda_{i} \lambda_{j} \gamma(\mathbf{s}_{i},\mathbf{s}_{j} ) +2 \sum\limits_{i=1}^{n}\lambda _{i} \gamma(\mathbf{s}_{i},\mathbf{s}_{0}) -2\sum\limits_{j=0}^{p}m_{j} \left( \sum\limits_{i=1}^{n}\lambda_{i} X_{j} (\mathbf{s}_{i})-X_{j} (\mathbf{s}_{0}) \right)\] Derivando respecto a \(\left\{ \lambda_{i} :i=1, \ldots,n\right\}\) y \(\left\{ m_{j} :j=0, \ldots,p\right\}\) e igualando a cero se obtienen las \(n+p+1\) ecuaciones del kriging universal que, expresadas en forma matricial, resultan ser: \[\boldsymbol{\Gamma}_{U} \boldsymbol{\lambda}_{U} =\boldsymbol{\gamma}_{U},\] con: \[\boldsymbol{\Gamma}_{U} = \left( \begin{array}{lc} \boldsymbol{\Gamma} & \mathbf{X} \\ \mathbf{X^\top } & \mathbf{0} \end{array} \right) ,\ \boldsymbol{\lambda}_{U} = \left( \begin{array}{c} \boldsymbol{\lambda} \\ \mathbf{m} \end{array} \right) ,\ \boldsymbol{\gamma}_{U} =\left( \begin{array}{c} \boldsymbol{\gamma} \\ \mathbf{x}_0 \end{array} \right),\] donde \(\boldsymbol{\gamma}=\left( \gamma(\mathbf{s}_{1},\mathbf{s}_{0}), \ldots, \gamma(\mathbf{s}_{n} ,\mathbf{s}_{0})\right)^\top\), \(\mathbf{m}=\left(m_{0}, \ldots,m_{p} \right)^\top\) y \(\boldsymbol{\Gamma}\) es una matriz \(n\times n\) con \(\boldsymbol{\Gamma}_{ij} = \gamma(\mathbf{s}_{i}, \mathbf{s}_{j})\). Además el MSPE mínimo, o varianza kriging: \[\sigma_{KU}^{2} (\mathbf{s}_{0})=2\sum\limits_{i=1}^{n}\lambda_{i} \gamma(\mathbf{s}_{0},\mathbf{s}_{i} )-\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}\lambda_{i} \lambda_{j} \gamma(\mathbf{s}_{i},\mathbf{s}_{j} )\] se puede obtener como: \[\begin{aligned} \sigma_{KU}^{2} (\mathbf{s}_{0}) & =\sum\limits_{i=1}^{n}\lambda_{i} \gamma(\mathbf{s}_{0},\mathbf{s}_{i})+\sum\limits_{j=0}^{p}m_{j} X_{j} (\mathbf{s}_{0}) \\ & =\boldsymbol{\lambda}_{U}^\top \boldsymbol{\gamma}_{U}. \end{aligned}\] En el caso particular del kriging ordinario (\(p=0\)), la expresión de la varianza kriging resulta ser: \[\sigma_{KO}^{2} (\mathbf{s}_{0})=\sum\limits_{i=1}^{n}\lambda_{i} \gamma(\mathbf{s}_{0},\mathbf{s}_{i})+m_{0}.\]

4.3.1 Ecuaciones en función del covariograma

Cuando existe el covariograma \(C(\mathbf{s}_{1},\mathbf{s}_{2}) = Cov(\varepsilon(\mathbf{s}_{1}), \varepsilon(\mathbf{s}_{2}))\) del proceso \(\varepsilon(\cdot)\) y es conocido (una suposición más fuerte), podemos expresar las ecuaciones del kriging universal (o del KO) en función de \(C(\cdot,\cdot)\). Además, si ninguna de las funciones explicativas es idénticamente 1, las ecuaciones del kriging universal sólo pueden expresarse en función del covariograma.

El proceso sería análogo al caso anterior, el sistema del kriging universal equivalente es: \[\boldsymbol{\Sigma}_{U} \boldsymbol{\lambda}_{U} = \mathbf{c}_{U},\] donde: \[\boldsymbol{\Sigma}_{U} =\left( \begin{array}{lc} \boldsymbol{\Sigma} & \mathbf{X} \\ \mathbf{X^\top } & \mathbf{0} \end{array} \right) ,\ \boldsymbol{\lambda}_{U} =\left( \begin{array}{c} \boldsymbol{\lambda} \\ \mathbf{m} \end{array} \right) ,\ \mathbf{c}_{U} =\left( \begin{array}{c} \mathbf{c} \\ \mathbf{x}_0 \end{array} \right),\] y la varianza kriging es: \[\begin{aligned} \sigma_{KU}^{2} (\mathbf{s}_{0}) & = C(\mathbf{s}_{0},\mathbf{s}_{0}) - \sum\limits_{i=1}^{n}\lambda_{i} C(\mathbf{s}_{0},\mathbf{s}_{i}) + \sum\limits_{j=0}^{p}m_{j} X_{j}(\mathbf{s}_{0}) \\ & = C(\mathbf{s}_{0}, \mathbf{s}_{0}) - \boldsymbol{\lambda}_{U} \mathbf{c}_{U}. \end{aligned}\]

Muchos de los algoritmos utilizados para la solución de los sistema kriging están diseñados y optimizados para covariogramas (e.g. Chilès y Delfiner, 2012, p. 170). En el caso de variogramas no acotados se podrían emplear pseudo-covarianzas (ver Sección 3.3.2).

4.3.2 Kriging residual

Otra forma, que puede ser más interesante, de obtener las ecuaciones del KU es a partir del predictor del kriging simple. Suponiendo que \(\boldsymbol{\beta}\) es conocido en el modelo del KU, el predictor del kriging simple es: \[\begin{aligned} p_{KS} (\mathbf{Z},\mathbf{s}_{0}) & = \mathbf{x}_0^{\top}\boldsymbol{\beta} + \mathbf{c}^{\top} \boldsymbol{\Sigma}^{-1} \left( \mathbf{Z} - \mathbf{X}\boldsymbol{\beta} \right) \\ & =\mathbf{c^\top }\boldsymbol{\Sigma}^{-1} \mathbf{Z}+(\mathbf{x}_0-\mathbf{X^\top }\boldsymbol{\Sigma}^{-1} \mathbf{c})^\top \boldsymbol{\beta}. \end{aligned}\]

Cuando \(\boldsymbol{\beta}\) no es conocido es lógico pensar en utilizar en su lugar su estimador lineal óptimo \[\hat{\boldsymbol{\beta}}_{gls} =(\mathbf{X}^{\top}\boldsymbol{\Sigma}^{-1} \mathbf{X})^{-1} \mathbf{X}^{\top}\boldsymbol{\Sigma}^{-1} \mathbf{Z},\] obteniéndose el predictor: \[p^{\ast}(\mathbf{Z},\mathbf{s}_{0}) = \mathbf{c^\top }\boldsymbol{\Sigma}^{-1}\mathbf{Z} + (\mathbf{x}_0-\mathbf{X}^\top\boldsymbol{\Sigma}^{-1} \mathbf{c})^\top \hat{\boldsymbol{\beta}}_{gls}.\] Puede verse (Goldberger, 1962) que este predictor, lineal e insesgado, es óptimo (en el sentido de que minimiza el MSPE sobre todos los predictores lineales e insesgados) y por tanto coincide con el predictor del kriging universal. Además, teniendo en cuenta que el error \(\left( p_{KS} (\mathbf{Z},\mathbf{s}_{0})-Z(\mathbf{s}_{0})\right)\) tiene covarianza nula con cualquier combinación lineal de \(\mathbf{Z}\) (ver e.g. Chilès y Delfiner, 2012, p. 161), esta relación también se extiende a la varianza kriging: \[\sigma_{KU}^{2} (\mathbf{s}_{0})=\sigma_{KS}^{2} (\mathbf{s}_{0} )+(\mathbf{x}_0-\mathbf{X^\top }\boldsymbol{\Sigma}^{-1} \mathbf{c})^\top \left( \mathbf{X^\top }\boldsymbol{\Sigma}^{-1} \mathbf{X}\right)^{-1} (\mathbf{x}_0-\mathbf{X^\top }\boldsymbol{\Sigma}^{-1} \mathbf{c}),\] donde el segundo termino cuantifica la precisión en la estimación de la media. Estas expresiones son conocidas como la relación de aditividad entre el KS y el KU.

Los resultados anteriores permiten pensar en la predicción lineal con media desconocida como un proceso de dos etapas: en la primera estimar la media desconocida, y en la segunda realizar la predicción lineal óptima con media supuestamente conocida. En el caso de una tendencia lineal obtenemos el predictor de KU, mientras que en el caso general se obtiene el denominado predictor del kriging residual (o kriging con tendencia externa).