# Kriging

Unter Kriging (oder auch: Krigen) versteht man ein geostatistisches Prognose- und Interpolationsverfahren, mit dem man eine räumlich verortete Variable an Orten, an denen sie nicht gemessen wurde, durch umliegende Messwerte interpolieren oder auch annähern kann. Stark vereinfacht könnte man sagen, diese Prognose ist eine Art gewichteter Mittelwert aus allen oder einigen der bekannten Messwerte einer Stichprobe. Außerhalb der Geostatistik ist das Verfahren als Gaußprozess-Regression bekannt.
Der südafrikanische Bergbauingenieur Danie Krige entwickelte 1951 für den Goldbergbau eine Interpolationsmethode, die auf der Abhängigkeit der Messwerte von den Abständen basiert, die zwischen den zugehörigen Messpunkten liegen. Der französische Mathematiker und Ingenieur Georges Matheron veröffentlichte 1960 die Arbeit „Krigeage d’un Panneau Rectangulaire par sa Périphérie“, welche die theoretische Grundlage der von Danie Krige entwickelten Methode schuf und sie nach ihm benannte.
Der wesentliche Vorteil gegenüber einfacheren Methoden wie beispielsweise der Inversen Distanzwichtung ist die Berücksichtigung der räumlichen Varianz, die sich mit Hilfe von Semivariogrammen ermitteln lässt. Die Semivarianz beschreibt, wie die Unterschiede zwischen den Messwerten zunehmen bzw. die Ähnlichkeit zwischen den Messwerten abnimmt, wenn der Abstand zwischen den Messpunkten größer wird. Sie eignet sich also dafür, die Gewichte der Mittelwertsbildung zu bestimmen, indem sie für näher gelegene Stichprobenwerte größere Gewichte, und für entferntere Stichprobenwerte kleinere Gewichte vergibt. Für einen gesuchten Wert werden dabei die Gewichte der in die Berechnung einfließenden Messwerte so bestimmt, dass der Prognosefehler möglichst gering ist. Der Prognosefehler hängt dabei von der Qualität des Variogramms bzw. der Variogrammfunktion ab, also wie gut das Semivariogrammmodell die tatsächliche räumliche Autokorrelation beschreibt.
Bei einfacheren Interpolationsverfahren können bei Häufung der Messpunkte Probleme auftreten. Dies wird beim Kriging vermieden und zwar durch die Berücksichtigung der statistischen Abstände zwischen der in die Berechnung eines Punktes einfließenden Nachbarn und Optimierung der gewichteten Mittel. Tritt an einer Stelle eine Clusterung auf, werden die Gewichte der Punkte innerhalb dieses Clusters gesenkt.
Unter Kriging versteht man die Bestimmung der besten linearen Prognose (oder Vorhersage) (englischer Kurzbegriff BLP, best linear prediction) eines Messwertes an einem nicht beobachteten Ort auf der Basis von Messwerten an beobachteten Orten. Dabei werden die Messwerte als Realisierungen von Zufallsvariablen (zufälliger Messungen) modelliert, die ein Zufallsfeld mit bekannter Erwartungswertfunktion und Kovarianzfunktion bilden. Werden in einem allgemeineren Kontext unbekannte Parameter in der Erwartungswertfunktion unverzerrt (erwartungstreu) geschätzt, so ergibt sich das Konzept der besten linearen unverzerrten Prognose (englischer Kurzbegriff BLUP, best linear unbiased prediction).

## Modellannahme

Die Messwerte $y_{0},y_{1},\dots ,y_{n}$ eines interessierenden Merkmals an den Orten $x_{0},x_{1},\dots ,x_{n}$ werden als Realisierungen von reellwertigen Zufallsvariablen $Y(x_{0}),Y(x_{1}),\dots ,Y(x_{n})$ aufgefasst, die ein Zufallsfeld $(Y(x))_{x\in {\mathcal {X}}}$ mit der Indexmenge ${\mathcal {X}}=\{x_{0},x_{1},\dots ,x_{n}\}$ bilden. Eigenschaften der gemeinsamen Wahrscheinlichkeitsverteilung des Zufallsfeldes sind die Erwartungswertfunktion
    $\mu (x)=\mathbb {E} [Y(x)],\quad x\in {\mathcal {X}}$
und die Kovarianzfunktion
    $k(x,z)=\mathrm {Cov} [(Y(x),Y(z)],\quad x,z\in {\mathcal {X}}\;.$
Im Spezialfall eines gaußschen Zufallsfeldes liegt durch die Erwartungswertfunktion und die Kovarianzfunktion die multivariate Wahrscheinlichkeitsverteilung der Zufallsvariablen $Y(x_{0}),Y(x_{1}),\dots ,Y(x_{n})$ als multivariate Normalverteilung fest.

## Statistische Fragestellung und Kriging-Lösung

Das Kriging beantwortet die Aufgabenstellung, einen Prognosewert ${\hat {y}}(x_{0})$ an einem nicht beobachteten Ort $x_{0}$ auf der Basis beobachteter Messwerte $y(x_{1}),\dots ,y(x_{n})$ an den Orten $x_{1},\dots ,x_{n}$ anzugeben. Im einfachsten Fall, wenn die Erwartungswertfunktion bekannt ist, handelt es sich um ein Prognoseverfahren. In Fällen, in denen Parameter der Ewartungswertfunktion geschätzt werden müssen, handelt es sich um eine kombiniertes Schätz- und Prognoseverfahren.

### Bekannte Parameter

Der einfachste Fall liegt vor, wenn die Erwartungswertfunktion und die Kovarianzfunktion bekannt sind. In diesem Fall ist die beste lineare Prognose (BLP) durch
    ${\hat {y}}(x_{0})=\mu _{0}+(\mathbf {y} -{\boldsymbol {\mu }})^{T}\mathbf {K} ^{-1}\mathbf {k}$
gegeben. Dabei gelten die folgenden Bezeichnungen:
* $\mu _{0}=\mathbb {E} [Y(x_{0})]$ bezeichnet den Erwartungswert der Zufallsvariablen $Y(x_{0})$,
* $\mathbf {y} =(y(x_{1}),\dots ,y(x_{n}))^{T}$ bezeichnet den Vektor der beobachteten Messwerte,
* ${\boldsymbol {\mu }}=(\mu _{1},\dots ,\mu _{n})^{T}$ ist ein ($n\times 1$)-Vektor mit den Erwartungswerten $\mu _{i}=\mathbb {E} [Y(x_{i})]$ für $i=1,\dots ,n$,
* $\mathbf {K} =[k_{ij}]_{i,j=1,\dots ,m}$ ist eine ($n\times n$)-Kovarianzmatrix mit den Kovarianzen $k_{ij}=\mathrm {Cov} [Y(x_{i}),Y(x_{j})]$ und
* $\mathbf {k} =(k_{1},\dots ,k_{n})^{T}$ ist ein ($n\times 1$)-Vektor mit den Kovarianzen $k_{i}=\mathrm {Cov} [Y(x_{0}),Y(x_{i})]$ für $i=1,\dots ,n$.
Der hochgestellte Index $T$ bezeichnet die Transponierung und $\mathbf {K} ^{-1}$ bezeichnet die invertierte Matrix.

### Geschätzte Parameter

Wenn die Erwartungswertfunktion unbekannt ist und durch eine bekannte Funktion $m:{\mathcal {X}}\to \mathbb {R} ^{p}$ und einen unbekannten Parametervektor ${\boldsymbol {\beta }}\in \mathbb {R} ^{p}$ in der Form
    $\mu (x)=m(x)^{T}{\boldsymbol {\beta }},\quad x\in {\mathcal {X}}$
dargestellt werden kann, ergibt sich ein kombiniertes Schätz- und Prognoseproblem, das durch die Angabe der besten linearen unverzerrten Prognose (BLUP) gelöst werden kann. Diese führt dann zum Prognosewert
    ${\hat {y}}(x_{0})=m(x_{0})^{T}{\hat {\boldsymbol {\beta }}}+(\mathbf {y} -\mathbf {M} {\boldsymbol {\hat {\beta }}})^{T}\mathbf {K} ^{-1}\mathbf {k}$
mit dem Schätzwert
    ${\hat {\boldsymbol {\beta }}}=(\mathbf {M} ^{T}\mathbf {K} ^{-1}\mathbf {M} )^{-1}\mathbf {M} ^{T}\mathbf {K} ^{-1}\mathbf {y}$
für den unbekannten Parametervektor ${\boldsymbol {\beta }}$. Dabei ist $\mathbf {M} :=[m(x_{1})\cdots m(x_{n})]^{T}$ eine ($n\times p$)-Matrix, bei der in der $i$-ten Zeile der Vektor $m(x_{i})^{T}$ steht.

## Beste lineare Prognose (BLP)
### Minimierung des mittleren quadratischen Prognosefehlers

Die beste lineare Prognose beruht auf der Minimierung des mittleren quadratischen Prognosefehlers
    $\mathbb {E} \left[\left({\hat {Y}}(x_{0})-Y(x_{0})\right)^{2}\right]$
zwischen der Prognose ${\hat {Y}}(x_{0})$ und der nicht beobachteten Zufallsvariable $Y(x_{0})$. Dabei wird bei der Minimierung die Menge der zulässigen Funktionen auf lineare Funktionen der Form
    ${\hat {Y}}(x_{0})=a_{0}+\sum _{i=1}^{n}a_{i}Y(x_{i}),\quad a_{0},a_{1},\dots ,a_{n}\in \mathbb {R}$
beschränkt. Die $Y(x_{i})$ sind die zufälligen Messungen an den Stellen $x_{1},\dotsc ,x_{n}$, $a_{0}$ ist eine Konstante ist und die Koeffizienten $a_{1},\dots ,a_{n}$ sind Gewichte der einzelnen Messungen. Die Bildung des Erwartungswertes bezieht sich sowohl auf die Zufallsvariable $Y(x_{0})$ als auch auf die Zufallsvariablen $Y(x_{1}),\dots ,Y(x_{n})$, die in ${\hat {Y}}(x_{0})$ eingehen.
Die Parameter $\lambda _{0},\lambda _{1},\dots ,\lambda _{n}$ der besten linearen Prognose im Sinn der Minimierung des mittleren quadratischen Prognosefehlers sind die Komponenten der Minimalstelle $(\lambda _{0},\lambda _{1},\dots ,\lambda _{n})$ mit der Eigenschaft
    $\mathbb {E} \left[\left(\lambda _{0}+\sum _{i=1}^{n}\lambda _{i}Y(x_{i})-Y(x_{0})\right)^{2}\right]=\min _{(a_{0},a_{1},\dots ,a_{n})^{T}\in \mathbb {R} ^{n+1}}\mathbb {E} \left[\left(a_{0}+\sum _{i=1}^{n}a_{i}Y(x_{i})-Y(x_{0})\right)^{2}\right]\;.$
Sie hängen nur von der Erwartungswertfunktion und der Kovarianzfunktion des Zufallsfeldes $(Y(x))_{x\in {\mathcal {X}}}$ ab. Für den Parametervektor ${\boldsymbol {\lambda }}=(\lambda _{1},\dots ,\lambda _{n})^{T}$ gilt
    ${\boldsymbol {\lambda }}=\mathbf {K} ^{-1}\mathbf {k}$
und für den Parameter $\lambda _{0}$ gilt
    $\lambda _{0}=\mu _{0}-{\boldsymbol {\mu }}^{T}{\boldsymbol {\lambda }}$.
Aus den bekannten Parametern $\mu _{0}$, ${\boldsymbol {\mu }}$, $\mathbf {k}$ und $\mathbf {K}$ können die Parameter $\lambda _{0}$ und ${\boldsymbol {\lambda }}$ bestimmt werden.
Der Minimalwert des mittleren quadratischen Prognosefehlers ist $k_{0}-\mathbf {k} ^{T}\mathbf {K} ^{-1}\mathbf {k}$ mit $k_{0}=\mathrm {Var} [Y(x_{0})]=k(Y(x_{0}),Y(x_{0}))$.

### Bester linearer Prognoswert

Für gegebene Messwerte $\mathbf {y} =(y(x_{1}),\dots ,y(x_{n}))^{T}$ als Realisierungen der Zufallsvariablen $Y(x_{1}),\dots ,Y(x_{n})$ ist dann der Wert
    ${\hat {y}}(x_{0})=\lambda _{0}+\sum _{i=1}^{n}\lambda _{i}y(x_{i})=\mu _{0}-{\boldsymbol {\mu }}^{T}{\boldsymbol {\lambda }}+\mathbf {y} ^{T}{\boldsymbol {\lambda }}=\mu _{0}+(\mathbf {y} -{\boldsymbol {\mu }})^{T}\mathbf {K} ^{-1}\mathbf {k}$
der beste lineare Prognosewert für die Messung am Ort $x_{0}$. Der beste lineare Prognosewert ist eine Realisierung der Zufallsvariablen
    ${\hat {Y}}(x_{0})=\lambda _{0}+\sum _{i=1}^{n}\lambda _{i}Y(x_{i})\;,$
die als beste lineare Prognose bezeichnet wird. Der zufällige Prognosefehler (die zufällige Prognoseabweichung) ${\hat {Y}}(x_{0})-Y(x_{0})$ hat den Erwartungswert 0 und die Varianz $k_{0}-\mathbf {k} ^{T}\mathbf {K} ^{-1}\mathbf {k}$.