UnterKriging (oder auch: Krigen) versteht man eingeostatistisches Prognose- und Interpolationsverfahren, mit dem man eine räumlich verorteteVariable an Orten, an denen sie nicht gemessen wurde, durch umliegendeMesswerteinterpolieren oder auchannähern kann. Stark vereinfacht könnte man sagen, diese Prognose ist eine Artgewichteter Mittelwert aus allen oder einigen der bekannten Messwerte einerStichprobe. Außerhalb der Geostatistik ist das Verfahren alsGaußprozess-Regression bekannt.^[1]

Der südafrikanische BergbauingenieurDanie 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 IngenieurGeorges Matheron veröffentlichte 1960 die Arbeit „Krigeage d’un Panneau Rectangulaire par sa Périphérie“^[2], 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 derInversen Distanzwichtung ist die Berücksichtigung der räumlichenVarianz, die sich mit Hilfe vonSemivariogrammen 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ächlicherä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 eineClusterung auf, werden die Gewichte der Punkte innerhalb dieses Clusters gesenkt.

Unter Kriging versteht man die Bestimmung derbesten linearen Prognose (oder Vorhersage) (englischer KurzbegriffBLP,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 einZufallsfeld mit bekannterErwartungswertfunktion undKovarianzfunktion bilden. Werden in einem allgemeineren Kontext unbekannte Parameter in der Erwartungswertfunktionunverzerrt (erwartungstreu) geschätzt, so ergibt sich das Konzept derbesten linearen unverzerrten Prognose (englischer KurzbegriffBLUP,best linear unbiased prediction).

Die Messwerte${\displaystyle y_0,y_1,\dots ,y_n}$ eines interessierenden Merkmals an den Orten${\displaystyle x_0,x_1,\dots ,x_n}$ werden als Realisierungen von reellwertigenZufallsvariablen${\displaystyle Y(x_0),Y(x_1),\dots ,Y(x_n)}$ aufgefasst, die einZufallsfeld${\displaystyle (Y(x){)}_{x\in \mathcal{X}}}$ mit der Indexmenge${\displaystyle \mathcal{X}=\{x_0,x_1,\dots ,x_n\}}$ bilden.^[3] Eigenschaften dergemeinsamen Wahrscheinlichkeitsverteilung des Zufallsfeldes sind dieErwartungswertfunktion

${\displaystyle \mu (x)=\mathbb{E}[Y(x)],x\in \mathcal{X}}$

und dieKovarianzfunktion

${\displaystyle k(x,z)=\mathrm{C}\mathrm{o}\mathrm{v}[(Y(x),Y(z)],x,z\in \mathcal{X}.}$

Im Spezialfall einesgaußschen Zufallsfeldes liegt durch die Erwartungswertfunktion und die Kovarianzfunktion diemultivariate Wahrscheinlichkeitsverteilung der Zufallsvariablen${\displaystyle Y(x_0),Y(x_1),\dots ,Y(x_n)}$ alsmultivariate Normalverteilung fest.

Das Kriging beantwortet die Aufgabenstellung, einen Prognosewert${\displaystyle \hat{y}(x_0)}$ an einemnicht beobachteten Ort${\displaystyle x_0}$ auf der Basis beobachteter Messwerte${\displaystyle y(x_1),\dots ,y(x_n)}$ an den Orten${\displaystyle 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.

Der einfachste Fall liegt vor, wenn die Erwartungswertfunktion und die Kovarianzfunktion bekannt sind. In diesem Fall ist diebeste lineare Prognose (BLP) durch

${\displaystyle \hat{y}(x_0)={\mu }_0+(\mathbf{y}-\mathit{\mu }{)}^T{\mathbf{K}}^{-1}\mathbf{k}}$

gegeben. Dabei gelten die folgenden Bezeichnungen:

• ${\displaystyle {\mu }_0=\mathbb{E}[Y(x_0)]}$ bezeichnet den Erwartungswert der Zufallsvariablen${\displaystyle Y(x_0)}$,
• ${\displaystyle \mathbf{y}=(y(x_1),\dots ,y(x_n){)}^T}$ bezeichnet den Vektor der beobachteten Messwerte,
• ${\displaystyle \mathit{\mu }=({\mu }_1,\dots ,{\mu }_n{)}^T}$ ist ein (${\displaystyle n\times 1}$)-Vektor mit den Erwartungswerten${\displaystyle {\mu }_i=\mathbb{E}[Y(x_i)]}$ für${\displaystyle i=1,\dots ,n}$,
• ${\displaystyle \mathbf{K}=[k_{ij}{]}_{i,j=1,\dots ,m}}$ ist eine (${\displaystyle n\times n}$)-Kovarianzmatrix mit den Kovarianzen${\displaystyle k_{ij}=\mathrm{C}\mathrm{o}\mathrm{v}[Y(x_i),Y(x_j)]}$ und
• ${\displaystyle \mathbf{k}=(k_1,\dots ,k_n{)}^T}$ ist ein (${\displaystyle n\times 1}$)-Vektor mit den Kovarianzen${\displaystyle k_i=\mathrm{C}\mathrm{o}\mathrm{v}[Y(x_0),Y(x_i)]}$ für${\displaystyle i=1,\dots ,n}$.

Der hochgestellte Index${\displaystyle T}$ bezeichnet dieTransponierung und${\displaystyle {\mathbf{K}}^{-1}}$ bezeichnet die invertierte Matrix.

Wenn die Erwartungswertfunktion unbekannt ist und durch eine bekannte Funktion${\displaystyle m:\mathcal{X}\to {\mathbb{R}}^p}$ und einen unbekannten Parametervektor${\displaystyle \mathit{\beta }\in {\mathbb{R}}^p}$ in der Form

${\displaystyle \mu (x)=m(x{)}^T\mathit{\beta },x\in \mathcal{X}}$

dargestellt werden kann, ergibt sich ein kombiniertes Schätz- und Prognoseproblem, das durch die Angabe derbesten linearen unverzerrten Prognose (BLUP) gelöst werden kann. Diese führt dann zum Prognosewert

${\displaystyle \hat{y}(x_0)=m(x_0{)}^T\hat{\mathit{\beta }}+(\mathbf{y}-\mathbf{M}\hat{\mathit{\beta }}{)}^T{\mathbf{K}}^{-1}\mathbf{k}}$

mit dem Schätzwert

${\displaystyle \hat{\mathit{\beta }}=({\mathbf{M}}^T{\mathbf{K}}^{-1}\mathbf{M}{)}^{-1}{\mathbf{M}}^T{\mathbf{K}}^{-1}\mathbf{y}}$

für den unbekannten Parametervektor${\displaystyle \mathit{\beta }}$. Dabei ist${\displaystyle \mathbf{M}:=[m(x_1)\cdots m(x_n){]}^T}$ eine (${\displaystyle n\times p}$)-Matrix, bei der in der${\displaystyle i}$-ten Zeile der Vektor${\displaystyle m(x_i{)}^T}$ steht.

Die beste lineare Prognose beruht auf der Minimierung desmittleren quadratischen Prognosefehlers

${\displaystyle \mathbb{E}\left[{\left(\hat{Y}(x_0)-Y(x_0)\right)}^2\right]}$

zwischen der Prognose${\displaystyle \hat{Y}(x_0)}$ und der nicht beobachteten Zufallsvariable${\displaystyle Y(x_0)}$.^[3] Dabei wird bei der Minimierung die Menge der zulässigen Funktionen auf lineare Funktionen der Form

${\displaystyle \hat{Y}(x_0)=a_0+\sum _{i=1}^n{a}_{i}Y(x_i),a_0,a_1,\dots ,a_n\in \mathbb{R}}$

beschränkt.^[3] Die${\displaystyle Y(x_i)}$ sind die zufälligen Messungen an den Stellen${\displaystyle x_1,\dots ,x_n}$,${\displaystyle a_0}$ ist eine Konstante ist und die Koeffizienten${\displaystyle a_1,\dots ,a_n}$ sind Gewichte der einzelnen Messungen. Die Bildung des Erwartungswertes bezieht sich sowohl auf die Zufallsvariable${\displaystyle Y(x_0)}$ als auch auf die Zufallsvariablen${\displaystyle Y(x_1),\dots ,Y(x_n)}$, die in${\displaystyle \hat{Y}(x_0)}$ eingehen.

Die Parameter${\displaystyle {\lambda }_0,{\lambda }_1,\dots ,{\lambda }_n}$ der besten linearen Prognose im Sinn der Minimierung des mittleren quadratischen Prognosefehlers sind die Komponenten der Minimalstelle${\displaystyle ({\lambda }_0,{\lambda }_1,\dots ,{\lambda }_n)}$ mit der Eigenschaft

${\displaystyle \mathbb{E}\left[{\left({\lambda }_0+\sum _{i=1}^n{\lambda }_{i}Y(x_i)-Y(x_0)\right)}^2\right]=\underset{(a_0,a_1,\dots ,a_n{)}^T\in {\mathbb{R}}^{n+1}}{min}\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${\displaystyle (Y(x){)}_{x\in \mathcal{X}}}$ ab.^[3] Für den Parametervektor${\displaystyle \mathit{\lambda }=({\lambda }_1,\dots ,{\lambda }_n{)}^T}$ gilt^[4]

${\displaystyle \mathit{\lambda }={\mathbf{K}}^{-1}\mathbf{k}}$

und für den Parameter${\displaystyle {\lambda }_0}$ gilt^[4]

${\displaystyle {\lambda }_0={\mu }_0-{\mathit{\mu }}^T\mathit{\lambda }}$.

Aus den bekannten Parametern${\displaystyle {\mu }_0}$,${\displaystyle \mathit{\mu }}$,${\displaystyle \mathbf{k}}$ und${\displaystyle \mathbf{K}}$ können die Parameter${\displaystyle {\lambda }_0}$ und${\displaystyle \mathit{\lambda }}$ bestimmt werden.

Der Minimalwert des mittleren quadratischen Prognosefehlers ist${\displaystyle k_0-{\mathbf{k}}^T{\mathbf{K}}^{-1}\mathbf{k}}$ mit${\displaystyle k_0=\mathrm{V}\mathrm{a}\mathrm{r}[Y(x_0)]=k(Y(x_0),Y(x_0))}$.^[4]

Für gegebene Messwerte${\displaystyle \mathbf{y}=(y(x_1),\dots ,y(x_n){)}^T}$ als Realisierungen der Zufallsvariablen${\displaystyle Y(x_1),\dots ,Y(x_n)}$ ist dann der Wert

${\displaystyle \hat{y}(x_0)={\lambda }_0+\sum _{i=1}^n{\lambda }_{i}y(x_i)={\mu }_0-{\mathit{\mu }}^T\mathit{\lambda }+{\mathbf{y}}^T\mathit{\lambda }={\mu }_0+(\mathbf{y}-\mathit{\mu }{)}^T{\mathbf{K}}^{-1}\mathbf{k}}$

der beste lineare Prognosewert für die Messung am Ort${\displaystyle x_0}$.Der beste lineare Prognosewert ist eine Realisierung der Zufallsvariablen

${\displaystyle \hat{Y}(x_0)={\lambda }_0+\sum _{i=1}^n{\lambda }_{i}Y(x_i),}$

die als beste lineare Prognose bezeichnet wird.Derzufällige Prognosefehler (die zufällige Prognoseabweichung)${\displaystyle \hat{Y}(x_0)-Y(x_0)}$ hat den Erwartungswert 0 und die Varianz${\displaystyle k_0-{\mathbf{k}}^T{\mathbf{K}}^{-1}\mathbf{k}}$.

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

${\displaystyle \underset{(a_0,a_1,\dots ,a_n{)}^T\in {\mathbb{R}}^{n+1}}{min}\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${\displaystyle ({\lambda }_0^{(1)},{\lambda }_1^{(1)},\dots ,{\lambda }_n^{(1)})}$ und${\displaystyle ({\lambda }_0^{(2)},{\lambda }_1^{(2)},\dots ,{\lambda }_n^{(2)})}$ gilt in diesem Fall

${\displaystyle \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

${\displaystyle \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.^[4]

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

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

und

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

Wenn für das Zufallsfeld${\displaystyle (Y(x){)}_{x\in \mathcal{X}}}$ die Erwartungswertfunktionund die Kovarianzfunktion bekannt sind, kann der beste lineare Prognosewert (BLP) für den Ort${\displaystyle x_0}$ basierend auf Messwerten${\displaystyle y_1,\dots ,y_n}$ an den Orten${\displaystyle 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änkte 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 dannbeste lineare unverzerrte Prognose (BLUP).

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

• Die Erwartungswertfunktion ist konstant, d. h.

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

wobei der Parameter${\displaystyle \mu }$ unbekannt ist.

• Die Kovarianzfunktion ist

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

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

DieKorrelationsfunktion ist in diesem Fall durch

${\displaystyle \varrho (x,y)=e^{-d(x,z)}>0\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.

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

${\displaystyle \mathbb{E}[Y(x_i)]=\mu ,\mathrm{V}\mathrm{a}\mathrm{r}[Y(x_i)]={\sigma }^2\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${\displaystyle (Y(x_i),Y(x_j))}$ durch${\displaystyle e^{-d(x_i,x_j)}}$ gegeben ist.

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

${\displaystyle \overline{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${\displaystyle \mu }$. Dabei bezeichnet${\displaystyle {\mathbf{1}}_n=(1,\dots ,1{)}^T\in {\mathbb{R}}^n}$ den Einsvektor der Dimension${\displaystyle n}$. Dagegen ergibt sich im hier vorliegenden Fall korrelierter Beobachtungen mit der verallgemeinerte Methode der kleinsten Quadrate der Schätzwert

${\displaystyle \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${\displaystyle \mu }$, der im Allgemeinen nicht mit dem arithmetischen Mittelwert${\displaystyle \overline{y}}$ übereinstimmt. Das sich die Komponenten des Gewichtsvektors${\displaystyle ({\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${\displaystyle {\mathbf{1}}_n}$ unmittelbar verifiziert, handelt es sich um einengewogenen arithmetischen Mittelwert der beobachteten Werte${\displaystyle y(x_1),\dots ,y(x_n)}$.

Im Spezialfall eines gaußschen Zufallsfeldes ist diemultivariate Wahrscheinlichkeitsverteilung der Zufallsvariablen${\displaystyle Y(x_1),\dots ,Y(x_n)}$ durch die Erwartungswertfunktion und die Kovarianzfunktion alsmultivariate Normalverteilung vollständig festgelegt, so dass es möglich ist, den Parameter${\displaystyle \mu }$ bei gegebenen Werten${\displaystyle y(x_1),\dots ,y(x_n)}$ durch dieMaximum-Likelihood-Methode zu bestimmen.^[7]

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

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

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

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

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

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

${\displaystyle \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

${\displaystyle 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${\displaystyle \mu }$ zurückzuführen ist.^[8][9]

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${\displaystyle \mathbf{m}:\mathcal{X}\to {\mathbb{R}}^p}$ mit${\displaystyle p\in \mathbb{N}}$ und einen unbekannten Parametervektor${\displaystyle \mathit{\beta }\in {\mathbb{R}}^p}$ in der Form

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

dargestellt werden kann. Es ergibt sich dann aus dem Ziel, eine Prognose für${\displaystyle Y(x_0)}$ abzugeben, ein kombiniertes Schätz- und Prognoseproblem, da für eine Prognose implizit der unbekannte Parametervektor${\displaystyle \mathit{\beta }}$ geschätzt werden muss. Die beste lineare unverzerrte Prognose kann mit einem zweistufigen Vorgehen gewonnen werden, bei dem zunächst der Parametervektor${\displaystyle \mathit{\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.^[8]

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

${\displaystyle Y(x)=\mathbf{m}(x{)}^T\mathit{\beta }+\epsilon (x),x\in \mathcal{X}}$

gesehen, wobei das Zufallsfeld${\displaystyle (\epsilon (x){)}_{x\in \mathcal{X}}}$ die konstante Erwartungswertfunktion${\displaystyle \mathbb{E}[\epsilon (x)]=0}$ für alle${\displaystyle x\in \mathcal{X}}$ und dieselbe Kovarianzfunktion wie das Zufallsfeld${\displaystyle (Y(x){)}_{x\in \mathcal{X}}}$ hat. Hierbei handelt es sich zunächst nur um eine andere, inhaltlich äquivalente Schreibweise für das Zufallsfeld${\displaystyle (Y(x){)}_{x\in \mathcal{X}}}$, indem man${\displaystyle \epsilon (x):=Y(x)-\mu (x)}$ für alle${\displaystyle x\in \mathcal{X}}$ definiert. Diese Schreibweise macht es aber möglich, die Theorie linearer Regressionsmodelle mit korrelierten Fehlertermen anzuwenden und den Schätzwert

${\displaystyle \hat{\mathit{\beta }}=({\mathbf{M}}^T{\mathbf{K}}^{-1}\mathbf{M}{)}^{-1}{\mathbf{M}}^T{\mathbf{K}}^{-1}\mathbf{y}}$

für den unbekannten Parametervektor${\displaystyle \mathit{\beta }}$ als beste lineare unverzerrte Schätzung mit bekannter Kovarianzmatrix zu bestimmen. Dabei bezeichnet${\displaystyle \mathbf{M}:=[\mathbf{m}(x_1)\cdots \mathbf{m}(x_n){]}^T}$ eine (${\displaystyle n\times p}$)-Matrix, bei der in der${\displaystyle i}$-ten Zeile der Vektor${\displaystyle \mathbf{m}(x_i{)}^T}$ steht.^[8]

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

${\displaystyle \hat{y}(x_0)=\mathbf{m}(x_0{)}^T\hat{\mathit{\beta }}+(\mathbf{y}-\mathbf{M}\hat{\mathit{\beta }}{)}^T{\mathbf{K}}^{-1}\mathbf{k}.}$

Der mittlere quadratische Prognosefehler der besten linearen unverzerrten Prognose ist durch

${\displaystyle k_0-{\mathbf{k}}^T{\mathbf{K}}^{-1}\mathbf{k}+{\mathit{\gamma }}^T({\mathbf{M}}^T{\mathbf{K}}^{-1}\mathbf{M}{)}^{-1}\mathit{\gamma }}$

mit

${\displaystyle \mathit{\gamma }=\mathbf{m}(x_0)-{\mathbf{M}}^T{\mathbf{K}}^{-1}\mathbf{k}}$

gegeben, wobei der letzte Term auf die Schätzung des Parametervektors${\displaystyle \mathit{\beta }}$ zurückzuführen ist.^[8]

• Beimeinfachen Kriging (simple Kriging) ist die Erwartungswertfunktion konstant,

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

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

• Beimgewö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.
• Beimuniversalen Kriging (universal Kriging) ist die Erwartungswertfunktion nicht konstant und wird durch einenlinearen Regressionsansatz modelliert. In diesem Fall kommt die beste lineare unverzerrte Prognose zur Anwendung, wobeiRegressionsparameter mitgeschätzt werden.
• Unterbayesianischem Kriging (bayesian Kriging) versteht man ein Verfahren bei dem der Schritt der Parameterschätzung mit Hilfe bayesianischer Schätzverfahren durchgeführt wird.
• DasIndikator-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 derinversen Distanzgewichtung (oder Distanzwichtung) ist der Prognosewert dergewogene arithmetische Mittelwert

${\displaystyle \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

${\displaystyle g_i=\frac{1}{d(x_0,x_i)}\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

mit${\displaystyle k_0>0}$ und${\displaystyle \alpha >0}$ hat. In diesem Fall ist${\displaystyle \mathbf{k}}$ derNullvektor,${\displaystyle {\mathbf{K}}^{-1}}$ ist eine Diagonalmatrik mit den Diagonalelementen${\displaystyle 1/(\alpha d(x_0,x_i)}$ für${\displaystyle i=1,\dots ,n}$ und der Schätzwert${\displaystyle ({\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${\displaystyle \sum _{i=1}^n{g}_{i}y(x_i)/\sum _{i=1}^n{g}_i}$, da sich der Faktor${\displaystyle \alpha }$ herauskürzt.

Teilweise wird in anwendungsnahen Darstellungen zugunsten einer vereinfachten Terminologie das Interpolations- und Prognoseproblem mit Begriffen aus der statistischen Schätztheorie beschrieben.^[10] Dadurch verschwimmt der Unterschied zwischen dem Schätzen eines unbekannten Parameters durch eineSchä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.^[1][11][12] 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.^[13]

• 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. 

Commons: Kriging – Sammlung von Bildern, Videos und Audiodateien

1. ↑^abMohamed A. Bouhlel, Joaquim R. R. A. Martins:Gradient-enhanced kriging for high-dimensional problems. In:Engineering with Computers.Band 35,Nr. 1, 1. Januar 2019,ISSN 1435-5663,doi:10.1007/s00366-018-0590-x.  Siehe Abschnitt2.1 Conventional kriging
2. ↑Centre de Géosciences/Géostatistique, Publications & documentation. Abgerufen am 18. April 2024. 
3. ↑^abcdMichael L. Stein:Interpolation of Spatial Data – Some Theory for Kriging.Abschnitt 1.2Best linear Prediction, S. 2. 
4. ↑^abcdefMichael L. Stein:Interpolation of Spatial Data – Some Theory for Kriging.Abschnitt 1.2Best linear Prediction, S. 3. 
5. ↑Michael L. Stein:Interpolation of Spatial Data – Some Theory for Kriging.Abschnitt 1.4An example of a poor BLP, S. 6–9. 
6. ↑Michael L. Stein:Interpolation of Spatial Data – Some Theory for Kriging.Abschnitt 1.5Best linear unbiased prediction, S. 7–9. 
7. ↑Michael L. Stein:Interpolation of Spatial Data – Some Theory for Kriging.Abschnitt 6.4Likelihood Methods, S. 169–175. 
8. ↑^abcdeMichael L. Stein:Interpolation of Spatial Data – Some Theory for Kriging.Abschnitt 1.5Best linear unbiased prediction, S. 8. 
9. ↑Die im Beispiel angegebenen Formeln ergeben sich aus der allgemeineren Darstellung in M. L. Stein,Interpolation of Spatial Data – Some Theory for Kriging, S. 7–8 mit den folgenden Spezialisierungen der dort verwendeten Notation:${\displaystyle p=1}$,${\displaystyle \mathit{\beta }=\mu }$,${\displaystyle \mathbf{m}(\mathbf{x})=1}$,${\displaystyle \mathbf{M}={\mathbf{1}}_n}$.
10. ↑Jörg Benndorf:Angewandte Geodatenanalyse und -Modellierung – Eine Einführung in die Geostatistik für Geowissenschaftler und Geoingenieure. Springer Vieweg, Wiesbaden 2023,ISBN 978-3-658-39980-1,Kap. 7Geostatistische Verfahren zur räumlichen Interpolation - Kriging, S. 157–201,doi:10.1007/978-3-658-39981-8. 
11. ↑Carl Edward Rasmussen, Christopher K. I. Williams:Gaussian Processes for Machine Learning. MIT Press, Cambridge / London 2006,ISBN 0-262-18253-X,S. 30 (gaussianprocess.org [PDF]). 
12. ↑Robert B. Gramacy:Surrogates – Gaussian Process Modeling, Design, and Optimization for the Applied Siences (= Texts in Statistical Science). CRC Press, Boca Raton / London / New York 2020,ISBN 978-1-03-224255-2,S. 143 (gramacy.com [PDF]). 
13. ↑Carl Edward Rasmussen, Christopher K. I. Williams:Gaussian Processes for Machine Learning. MIT Press, Cambridge / London 2006,ISBN 0-262-18253-X,S. 16 (gaussianprocess.org [PDF]). 

• Geostatistik
• Ausgleichsrechnung