Kriging

Geländemodell vor und nach Kriging-Interpolation

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 ein 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}$.

Anmerkungen

• Die obige Darstellung enthält die stillschweigend gemachte Voraussetzung, dass die Zufallsvariablen $Y(x_{i})$ für $i=0,1,\dots ,n$ jeweils eine endliche Varianz besitzen. Ohne diese Annahme, die beispielsweise für Cauchy-verteilte Zufallsvariablen verletzt ist, ist das Kriterium der Minimierung des mittleren quadratischen Prognosefehlers nicht anwendbar.
• Es ist oben vorausgesetzt, dass die Kovarianzmatrix $\mathbf {K}$ invertierbar ist und somit die Kovarianzmatrix nicht nur positiv semidefinit, sondern positiv definit ist. Falls die Kovarianzmatrix nicht invertierbar ist, ist für

$\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]$

zwar das Minimum, nicht aber die Minimalstelle eindeutig. Für zwei verschiedene minimierende Parametervektoren $(\lambda _{0}^{(1)},\lambda _{1}^{(1)},\dots ,\lambda _{n}^{(1)})$ und$(\lambda _{0}^{(2)},\lambda _{1}^{(2)},\dots ,\lambda _{n}^{(2)})$ gilt in diesem Fall

$\mathbb {E} \left[\left(\lambda _{0}^{(1)}+\sum _{i=1}^{n}\lambda _{i}^{(1)}Y(x_{i})-\left(\lambda _{0}^{(2)}+\sum _{i=1}^{n}\lambda _{i}^{(2)}Y(x_{i}))\right)\right)^{2}\right]=0,$

sodass zwar die Parameter nicht eindeutig sind, aber der prognostizierte Wert

${\hat {y}}(x_{0})=\lambda _{0}^{(1)}+\sum _{i=1}^{n}\lambda _{i}^{(1)}y(x_{i})=\lambda _{0}^{(2)}+\sum _{i=1}^{n}\lambda _{i}^{(2)}y(x_{i})$

eindeutig ist.

• Die beste lineare Prognose ist im Allgemeinen nicht die beste Prognose im Sinn der Minimierung des mittleren quadratischen Prognosefehlers, wenn für die Prognose ${\hat {Y}}(x_{0})$ allgemeinere nichtlineare Funktionen $g(Y(x_{1}),\dots ,Y(x_{n}))$ zugelassen werden.
• Eine Besonderheit ergibt sich im Spezialfall eines gaußschen Zufallsfelds. Bei dieser speziellen Verteilungsannahme ist die beste Prognose bezüglicher aller, auch nicht-linearer Funktionen, durch die beste lineare Prognose gegeben. In diesem Fall ist also die beste lineare Prognose zugleich die beste Prognose. Außerdem ist die bedingte Verteilung von $Y(x_{0})$ gegeben $(Y(x_{1}),\dots ,Y(x_{n}))=(y(x_{1}),\dots ,y(x_{n}))$ eine univariate Normalverteilung ${\mathcal {N}}(\mu _{0|1,\dots ,n},\sigma _{0|1,\dots ,n}^{2})$ mit den Parametern

$\mu _{0|1,\dots ,n}=\lambda _{0}+\sum _{i=1}^{n}\lambda _{i}y(x_{i})$

und

$\sigma _{0|1,\dots ,n}^{2}=k_{0}-\mathbf {k} ^{T}\mathbf {K} ^{-1}\mathbf {k} \;.$

Beste lineare unverzerrte Prognose (BLUP)

Wenn für das Zufallsfeld $(Y(x))_{x\in {\mathcal {X}}}$ die Erwartungswertfunktion und die Kovarianzfunktion bekannt sind, kann der beste lineare Prognosewert (BLP) für den Ort $x_{0}$ basierend auf Messwerten $y_{1},\dots ,y_{n}$ an den Orten $x_{1},\dots ,x_{n}$ einfach, wie im vorausgegangenen Abschnitt angegeben, berechnet werden.

Wenn aber, was eher typisch ist, die Erwartungswert- und Kovarianzfunktion teilweise unbekannt sind, sind verschiedene einschränkende Modellannahmen erforderlich, um die Zahl der unbekannten Parameter so zu senken, dass diese mit den vorhandenen Beobachtungen schätzbar sind. Es entsteht dann ein kombiniertes Schätz- und Prognoseproblem. Ein bestimmtes Verfahren, bei dem die unbekannten Parameter im Rahmen eines linearen Modellansatzes unverzerrt (erwartungstreu) geschätzt werden, heißt dann beste lineare unverzerrte Prognose (BLUP).

Beispiel

Ein einfaches Beispielmodell beruht auf den beiden folgenden stark vereinfachenden Annahmen:

• Die Erwartungswertfunktion ist konstant, d. h.

$\mu (x)=\mu \quad {\text{für alle }}x\in {\mathcal {X}}\;,$

wobei der Parameter $\mu$ unbekannt ist.

• Die Kovarianzfunktion ist

$k(x,z)=\sigma ^{2}e^{-d(x,z)}\quad {\text{für alle }}x,z\in {\mathcal {X}}\;,$

wobei der Parameter $\sigma ^{2}>0$ bekannt ist. Dabei bezeichnet $d(x,y)=\|x-z\|$ die euklidische Distanz zwischen den Orten $x$ und $z$. Die Distanz kann zweidimensional in der Fläche, im dreidimensionalen Raum oder allgemeiner in einem $d$-dimensionalen Raum gemessen werden, in dem die euklidische Distanz definiert ist.

Die Korrelationsfunktion ist in diesem Fall durch

$\varrho (x,y)=e^{-d(x,z)}>0\quad {\text{für alle }}x,z\in {\mathcal {X}}$

gegeben. Mit zunehmender Distanz nimmt die Korrelation ab. Mit abnehmender Distanz nähert sich die Korrelation dem Wert Eins. Da die Koordinaten der Orte als bekannt vorausgesetzt sind, können die Distanzen und damit die Werte der Kovarianzfunktion bestimmt werden.

Parameterschätzung

Im Beispielmodell gibt es den unbekannten Parameter $\mu$. Da die beiden Parameter $\mu$ und $\sigma ^{2}$ zugleich der Erwartungswert und die Varianz der $n$ beobachtbaren Zufallsvariablen $Y(x_{1}),\dots ,Y(x_{n})$ sind, also

$\mathbb {E} [Y(x_{i})]=\mu ,\quad \mathrm {Var} [Y(x_{i})]=\sigma ^{2}\quad {\text{für }}i=1,\dots ,n$

gilt, scheint bei oberflächlicher Betrachtung ein Standardproblem der statistischen Schätztheorie vorzuliegen. Dies ist aber nicht der Fall, da bei Standardproblemen der statistischen Schätztheorie von stochastisch unabhängigen Beobachtungen ausgegangen wird. In diesem Beispiel sind aber alle Beobachtungspaare positiv korreliert, wobei die Korrelation für ein Paar $(Y(x_{i}),Y(x_{j}))$ durch $e^{-d(x_{i},x_{j})}$ gegeben ist.

Die Schätzung des Parameters $\mu$ aus gegebenen beobachteten Werten $y(x_{1}),\dots ,y(x_{n})$ kann in diesem Modellzusammenhang mit Hilfe der verallgemeinerten Methode der kleinsten Quadrate erfolgen, die es ermöglicht, bei der Schätzung eine gegebene Korrelationsstruktur zwischen den beobachtbaren Variablen zu berücksichtigen. Im Fall unkorrelierter Variablen führt die gewöhnliche Methode der kleinsten Quadrate zu dem üblichen Schätzwert

${\bar {y}}={\frac {\mathbf {1} _{n}^{T}\mathbf {y} }{\mathbf {1} _{n}^{T}\mathbf {1} _{n}}}={\frac {1}{n}}\sum _{i=1}^{n}y_{i}$

für den Parameter $\mu$. Dabei bezeichnet $\mathbf {1} _{n}=(1,\dots ,1)^{T}\in \mathbb {R} ^{n}$ den Einsvektor der Dimension $n$. Dagegen ergibt sich im hier vorliegenden Fall korrelierter Beobachtungen mit der verallgemeinerte Methode der kleinsten Quadrate der Schätzwert

${\hat {\mu }}={\frac {\mathbf {1} _{n}^{T}\mathbf {K} ^{-1}\mathbf {y} }{\mathbf {1} _{n}^{T}\mathbf {K} ^{-1}\mathbf {1} _{n}}}=(\mathbf {1} _{n}^{T}\mathbf {K} ^{-1}\mathbf {1} _{n})^{-1}\mathbf {1} _{n}^{T}\mathbf {K} ^{-1}\mathbf {y}$

für den Parameter $\mu$, der im Allgemeinen nicht mit dem arithmetischen Mittelwert ${\bar {y}}$ übereinstimmt. Das sich die Komponenten des Gewichtsvektors $(\mathbf {1} _{n}^{T}\mathbf {K} ^{-1}\mathbf {1} _{n})^{-1}\mathbf {1} _{n}^{T}\mathbf {K} ^{-1}$ zu Eins addieren, wie man durch Rechtsmultiplikation mit dem Vektor $\mathbf {1} _{n}$ unmittelbar verifiziert, handelt es sich um einen gewogenen arithmetischen Mittelwert der beobachteten Werte $y(x_{1}),\dots ,y(x_{n})$.

Im Spezialfall eines gaußschen Zufallsfeldes ist die multivariate Wahrscheinlichkeitsverteilung der Zufallsvariablen $Y(x_{1}),\dots ,Y(x_{n})$ durch die Erwartungswertfunktion und die Kovarianzfunktion als multivariate Normalverteilung vollständig festgelegt, so dass es möglich ist, den Parameter $\mu$ bei gegebenen Werten $y(x_{1}),\dots ,y(x_{n})$ durch die Maximum-Likelihood-Methode zu bestimmen.

Prognose

Wäre der Parameter $\mu$ bekannt, so ergäbe sich der beste lineare Prognosewert für den Messwert $y(x_{0})$, wie oben angegeben, als

$\mu +(\mathbf {y} -\mu \mathbf {1} _{n})^{T}\mathbf {K} ^{-1}\mathbf {k} ,$

wobei der mittlere quadratische Prognosefehler durch $k_{0}-\mathbf {k} ^{T}\mathbf {K} ^{-1}\mathbf {k}$ gegeben ist.

Wenn in einem ersten Schritt ein Schätzwert ${\hat {\mu }}$ für den Parameter $\mu$ bestimmt ist, kann in einem zweiten Schritt mit der geschätzten Erwartungswertfunktion

${\hat {m}}(x)={\hat {\mu }}\quad {\text{für alle }}x\in {\mathcal {X}}\;,$

der beste lineare Prognosewert so bestimmt werden, als ob die Erwartungswertfunktion bekannt sei. Es ergibt sich dann

${\hat {y}}(x_{0})={\hat {\mu }}+(\mathbf {y} -{\hat {\mu }}\mathbf {1} _{n})^{T}\mathbf {K} ^{-1}\mathbf {k}$

als bester linearer unverzerrter Prognosewert. Der mittlere quadratische Prognosefehler erhöht sich durch den zufälligen Schätzfehler, der durch die Parameterschätzung verursacht ist. Der mittlere quadratische Prognosefehler der besten linearen unverzerrten Prognose ist durch

$k_{0}-\mathbf {k} ^{T}\mathbf {K} ^{-1}\mathbf {k} +{\frac {(1-\mathbf {1} _{n}^{T}\mathbf {K} ^{-1}\mathbf {k} )^{2}}{\mathbf {1} _{n}^{T}\mathbf {K} ^{-1}\mathbf {1} _{n}}}$

gegeben, wobei der letzte Term auf die Schätzung des Parameters $\mu$ zurückzuführen ist.

Allgemeiner Fall

Der allgemeine Fall des Kriging, in dem eine beste lineare unverzerrte Prognose durchführbar ist, liegt vor, wenn die Erwartungswertfunktion teilweise unbekannt ist und durch eine bekannte vektorwertige Funktion $\mathbf {m$ :{\mathcal {X}}\to \mathbb {R} ^{p}} mit $p\in \mathbb {N}$ und einen unbekannten Parametervektor ${\boldsymbol {\beta }}\in \mathbb {R} ^{p}$ in der Form

$\mu (x)=\mathbf {m} (x)^{T}{\boldsymbol {\beta }},\quad x\in {\mathcal {X}}$

dargestellt werden kann. Es ergibt sich dann aus dem Ziel, eine Prognose für $Y(x_{0})$ abzugeben, ein kombiniertes Schätz- und Prognoseproblem, da für eine Prognose implizit der unbekannte Parametervektor ${\boldsymbol {\beta }}$ geschätzt werden muss. Die beste lineare unverzerrte Prognose kann mit einem zweistufigen Vorgehen gewonnen werden, bei dem zunächst der Parametervektor ${\boldsymbol {\beta }}$ in einem linearen Modell geschätzt wird und dann mit den gewonnenen Schätzwerten formal eine beste lineare Prognose so berechnet wird, als ob die Schätzwerte die unbekannten Parameter seien.

Parameterschätzung

Der unbekannte Parametervektor ${\boldsymbol {\beta }}$ wird im Rahmen des linearen Modells

$Y(x)=\mathbf {m} (x)^{T}{\boldsymbol {\beta }}+\varepsilon (x),\quad x\in {\mathcal {X}}$

gesehen, wobei das Zufallsfeld $(\varepsilon (x))_{x\in {\mathcal {X}}}$ die konstante Erwartungswertfunktion $\mathbb {E} [\varepsilon (x)]=0$ für alle $x\in {\mathcal {X}}$ und dieselbe Kovarianzfunktion wie das Zufallsfeld $(Y(x))_{x\in {\mathcal {X}}}$ hat. Hierbei handelt es sich zunächst nur um eine andere, inhaltlich äquivalente Schreibweise für das Zufallsfeld $(Y(x))_{x\in {\mathcal {X}}}$, indem man $\varepsilon (x):=Y(x)-\mu (x)$ für alle $x\in {\mathcal {X}}$ definiert. Diese Schreibweise macht es aber möglich, die Theorie linearer Regressionsmodelle mit korrelierten Fehlertermen anzuwenden und den 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 }}$ als beste lineare unverzerrte Schätzung mit bekannter Kovarianzmatrix zu bestimmen. Dabei bezeichnet $\mathbf {M$ :=[\mathbf {m} (x_{1})\cdots \mathbf {m} (x_{n})]^{T}} eine ($n\times p$)-Matrix, bei der in der $i$-ten Zeile der Vektor $\mathbf {m} (x_{i})^{T}$ steht.

Prognose

Die besten linearen unverzerrten Prognose ergibt sich, wenn der Prognosewert analog zum Vorgehen bei der besten linearen Prognose mit bekannter Erwartungsfunktion so bestimmt wird, als ob der Schätzwert der Parameter wäre. Damit ergibt sich der Prognosewert

${\hat {y}}(x_{0})=\mathbf {m} (x_{0})^{T}{\hat {\boldsymbol {\beta }}}+(\mathbf {y} -\mathbf {M} {\hat {\boldsymbol {\beta }}})^{T}\mathbf {K} ^{-1}\mathbf {k} \;.$

Der mittlere quadratische Prognosefehler der besten linearen unverzerrten Prognose ist durch

$k_{0}-\mathbf {k} ^{T}\mathbf {K} ^{-1}\mathbf {k} +{\boldsymbol {\gamma }}^{T}(\mathbf {M} ^{T}\mathbf {K} ^{-1}\mathbf {M} )^{-1}{\boldsymbol {\gamma }}$

mit

${\boldsymbol {\gamma }}=\mathbf {m} (x_{0})-\mathbf {M} ^{T}\mathbf {K} ^{-1}\mathbf {k}$

gegeben, wobei der letzte Term auf die Schätzung des Parametervektors ${\boldsymbol {\beta }}$ zurückzuführen ist.

Spezialfälle des Kriging

• Beim einfachen Kriging (simple Kriging) ist die Erwartungswertfunktion konstant,

$\mu (x)=\mu$ für alle $x\in {\mathcal {X}}\;,$

und der Parameter $\mu$ ist bekannt. In diesem Fall kommt die beste lineare Prognose zur Anwendung.

• Beim gewöhnlichen Kriging (ordinary Kriging) ist die Erwartungswertfunktion konstant, aber der gemeinsame Erwartungswert ist unbekannt und muss aus den beobachteten Werten geschätzt werden. In diesem Fall kommt die beste lineare unverzerrte Prognose zur Anwendung. Das oben ausgeführte Beispiel ist ein Fall des gewöhnlichen Kriging.
• Beim universalen Kriging (universal Kriging) ist die Erwartungswertfunktion nicht konstant und wird durch einen linearen Regressionsansatz modelliert. In diesem Fall kommt die beste lineare unverzerrte Prognose zur Anwendung, wobei Regressionsparameter mitgeschätzt werden.
• Unter bayesianischem Kriging (bayesian Kriging) versteht man ein Verfahren bei dem der Schritt der Parameterschätzung mit Hilfe bayesianischer Schätzverfahren durchgeführt wird.
• Das Indikator-Kriging ist ein Spezialfall bei dem die beobachteten Werte nur die Werte 0 und 1 annehmen, beispielsweise den Wert 0, wenn ein Grenzwert nicht überschritten ist, und den Wert 1, wenn ein Grenzwert überschritten ist.
• Bei der inversen Distanzgewichtung (oder Distanzwichtung) ist der Prognosewert der gewogene arithmetische Mittelwert

${\hat {y}}(x_{0})={\frac {\sum _{i=1}^{n}g_{i}y(x_{i})}{\sum _{i=1}^{n}g_{i}}}$

der beobachteten Werte mit den positiven Gewichten

$g_{i}={\frac {1}{d(x_{0},x_{i})}}\quad {\text{für }}i=1,\dots ,n\;.$

Wie das obige Beispiel zeigt, ergibt sich dieser Fall als beste lineare unverzerrte Prognose, und damit als Spezialfall des gewöhnlichen Kriging, wenn die Erwartungswertfunktion konstant ist und die Kovarianzfunktion die spezielle Form

$k(x,z)={\begin{cases}k_{0}&{\text{für }}x=z=x_{0}\\\alpha \cdot d(x_{0},z)&{\text{für }}x=z\neq x_{0}\\0&{\text{für }}x\neq z\end{cases}},\quad x,z\in {\mathcal {X}}$

mit $k_{0}>0$ und $\alpha >0$ hat. In diesem Fall ist $\mathbf {k}$ der Nullvektor, $\mathbf {K} ^{-1}$ ist eine Diagonalmatrik mit den Diagonalelementen $1/(\alpha d(x_{0},x_{i})$ für $i=1,\dots ,n$ und der Schätzwert $(\mathbf {1} _{n}^{T}\mathbf {K} ^{-1}\mathbf {1} _{n})^{-1}\mathbf {1} _{n}^{T}\mathbf {K} ^{-1}\mathbf {y}$ vereinfacht sich zum gewogenen arithmetische Mittelwert $\sum _{i=1}^{n}g_{i}y(x_{i})/\sum _{i=1}^{n}g_{i}$, da sich der Faktor $\alpha$ herauskürzt.

Abweichende Interpretationen und Verallgemeinerungen

Teilweise wird in anwendungsnahen Darstellungen zugunsten einer vereinfachten Terminologie das Interpolations- und Prognoseproblem mit Begriffen aus der statistischen Schätztheorie beschrieben. Dadurch verschwimmt der Unterschied zwischen dem Schätzen eines unbekannten Parameters durch eine Schätzfunktion und der Prognose des Wertes einer Zufallsvariablen.

Im engeren Sinn bezeichnet Kriging die oben beschriebene Modellierungsmethode der zu schätzenden Parameter durch ein lineares Modell und die dann explizit angebbaren Lösungen im Sinn der besten linearen unverzerrten Prognose. Teilweise wird Kriging aber auch allgemeiner, orientiert an der Fragestellung des Kringing, für andere methodische Vorgehensweisen verwendet. So werden im Bereich des maschinellen Lernens Methoden der Gaußprozess-Regression basierend auf gaußschen Zufallsfeldern als Kriging bezeichnet. Die klassische Kriging-Methode benötigt keine Normalverteilungsannahme und verarbeitet nur die Informationen der Erwartungswert- und Kovarianzfunktion im Rahmen einer linearen Modellstruktur mit klassischen statistischen Methoden. Dagegen wird bei der Gaußprozess-Regression durch eine weitgehende Annahme einer multivariaten Normalverteilung für alle Zufallsvariablen die Möglichkeit eröffnet, für die multivariate Normalverteilung zur Verfügung stehende Methoden zur Bestimmung von Prognoseverteilungen als bedingten Wahrscheinlichkeitsverteilungen zu verwenden.

Literatur

• Danie G. Krige: A statistical approach to some basic mine valuation problems on the Witwatersrand. In: J. of the Chem., Metal. and Mining Soc. of South Africa. 52 (6), 1951, S. 119–139.
• Rudolf Dutter: Mathematische Methoden in der Technik. Band 2: Geostatistik. B.G. Teubner Verlag, Stuttgart 1985, ISBN 3-519-02614-7.
• J. P. Chiles, P. Delfiner: Geostatistics: Modeling Spatial Uncertainty. Wiley, New York 1999, ISBN 0-471-08315-1.
• Michael Leonhard Stein: Interpolation of Spatial Data – Some Theory for Kriging (= Springer Series in Statistics). Springer, New York 1999, ISBN 978-1-4612-7166-6, doi:10.1007/978-1-4612-1494-6. 

Weblinks

Commons: Kriging – Sammlung von Bildern, Videos und Audiodateien