Kriging as an Interpolation Method, or Gaussian Process Regression

In statistics, originally in geostatistics, kriging or Gaussian process regression is an interpolation method in which the interpolated values are modeled by a Gaussian process governed by prior covariances, as opposed to a piecewise-polynomial spline chosen to optimize the smoothness of the interpolated values. Under suitable assumptions about the priors, kriging gives the best linear unbiased prediction of the intermediate values. Interpolation methods based on other criteria, such as smoothness (e.g. a smoothing spline), may not yield the most likely intermediate values. The method is widely used in the field of spatial analysis and computer experiments. This method is also known as Wiener–Kolmogorov prediction, after Norbert Wiener and Andrey Kolmogorov.

Example of one-dimensional data interpolation by kriging, with confidence intervals. Squares indicate the location of the data. The kriging interpolation, shown in red, runs along the mean values of the normally distributed confidence intervals, shown in gray. The dashed curve shows a spline, which is smooth but deviates significantly from the expected intermediate values given by these mean values.

The theoretical foundation of the method was developed by the French mathematician Georges Matheron in 1960, based on the master's thesis of Danie G. Krige, a pioneer in constructing a plot of the weighted-average gold content over the Witwatersrand reef complex in South Africa. Krige sought to estimate the most likely distribution of gold based on samples from a few boreholes. The English verb is to krige, and the most common noun is kriging; both are frequently pronounced with a hard "g", following the English pronunciation of the name "Krige". In the literature, the word is sometimes capitalized.

Although in its basic formulation kriging requires substantial computational resources, it can be scaled to solve larger problems using various approximation methods.

This interpolation method is named after the South African mining engineer Daniel Krige (eng.)rus., who manually produced geological maps from a limited set of data over a given area. It is a form of generalized linear regression that uses statistical parameters to find the optimal estimate, in the sense of minimum mean-square deviation, when constructing surfaces, cubes, and maps. The method is based on the principle of unbiasedness of the mean; that is, taken together, the values on the map should have the correct mean value. Global unbiasedness is formally ensured by raising the low values and lowering the high ones.

Given the correctly chosen prior assumptions, kriging gives the best linear unbiased prediction of intermediate values. Interpolation methods based on other criteria, such as smoothness, need not give the most likely values at intermediate points. This method is widely used in the field of spatial analysis and computer (numerical) experiments. This method is also known as Wiener–Kolmogorov prediction, after Norbert Wiener and Andrey Nikolaevich Kolmogorov.

From the point of view of general statistics, kriging consists of minimizing the variance of the measurement error, which is a function of the measurement weights. Minimizing this variance reduces the mean square error of the deviation of the estimated value from the possible one. This is achieved by setting the first derivative of the error with respect to each unknown weight to zero. This ultimately yields a system of equations whose solution is the vector of weights.

Kriging performs two groups of tasks:

1. quantifying the spatial structure of the data
2. producing a prediction

The quantitative representation of the spatial structure of the data, known as variogram construction, allows users to fit a model of spatial dependence to the data. To calculate (predict) the unknown value of a variable at a given location, kriging uses the fitted variogram model, the configuration of the spatial data, and the values at the measurement points around the given location.

Basic principles

Related terms and methods

The basic idea of kriging is to predict the value of a function at a given point by computing a weighted average of the known values of the function in the neighborhood of the point. This method is mathematically closely related to regression analysis. Both theories derive the best linear unbiased estimator, based on assumptions about covariances, use the Gauss–Markov theorem to prove the independence of the estimate and the error, and use very similar formulas. Nevertheless, they are useful in different frameworks: kriging is used to estimate a single realization of a random field, whereas regression models are based on multiple observations of a multivariate data set.

The kriging estimate can also be viewed as a spline in a reproducing kernel Hilbert space, with the reproducing kernel given by the covariance function. The difference from the classical kriging approach lies in the interpretation: whereas the spline is motivated by minimum-norm interpolation based on the structure of the Hilbert space, kriging is motivated by the expected squared prediction error based on a stochastic model.

Kriging with polynomial trend surfaces is mathematically identical to generalized least-squares polynomial curve fitting.

Kriging can also be understood as a form of Bayesian inference. Kriging starts with a prior distribution over functions. This prior takes the form of a Gaussian process:samples from the function will be normally distributed, where the covariance between any two samples is the covariance function (or kernel) of the Gaussian process, evaluated at the spatial locations of the two points. A set of values is then observed, each value associated with a spatial location. A new value can now be predicted at any new spatial location by combining the Gaussian prior function with a Gaussian likelihood function for each of the observed values. The resulting posterior distribution is also Gaussian, with a mean and covariance that can simply be computed from the observed values, their variance, and the kernel matrix derived from the prior.

Geostatistical estimation

In geostatistical models, the sample data are interpreted as the result of a random process. The fact that these models incorporate uncertainty in their conceptualization does not mean that the phenomenon — a forest, an aquifer, a mineral deposit — arose as the result of a random process, but rather it allows a methodological framework to be built for spatial inference of a quantity at unobserved locations and for quantifying the uncertainty associated with the estimate.

In the context of this model, a stochastic process is simply a way of approaching a data set collected from samples. The first step in geostatistical modeling is to build the random process that best describes the observed data set.

The value at location (a general notation for a set of geographic coordinates) is interpreted as a realizationof the random variable . In the space, over which the set of samples is scattered, there are realizations of random variables , which are correlated with each other.

The set of random variables represents a random function, of which only a single realization is known. is the set of observed data. With only a single realization of each random variable, it is theoretically impossible to determine any statistical parameter of the individual variables or of the function. The solution proposed in the geostatistical formalism is to assume some degree of stationarity of the random function in order to make it possible to infer certain statistical values.

For example, if we assume, based on homogeneity of the samples over the area over which the variable is distributed, the hypothesis that the first moment is stationary (i.e. all random variables have the same mean value), then it is assumed that the mean can be estimated using the arithmetic mean of the sample values.

The stationarity hypothesis relating to the second moment is defined as follows: the correlation between two random variables depends solely on the spatial distance between them and not on their location. Thus, if and then:

and for simplicity we define and .

This hypothesis allows us to derive these two measures — the variogram and the covariogram:

where:

• ;
• denotes the set of observation pairs such that , and is the number of pairs in the set. In this set, and denote the same element. Usually an «approximate distance» is used, implemented using a certain tolerance.

Linear estimation

Spatial inference, or estimation of the quantity , at an unobserved location , is computed from a linear combination of the observed values and weights :

The weights are intended to summarize two extremely important procedures in the spatial inference process:

• they reflect the structural «closeness» of the samples to the estimation location,
• at the same time they must have a desegregation effect, to avoid the bias caused by possible sample clusters

When computing the weights , there are two objectives in the geostatistical formalism: unbiasedness and minimum estimation variance.

If the cloud of real values is plotted against the estimated values , the criterion of global unbiasedness, intrinsic stationarity, or wide-sense stationarity of the field, implies that the mean of the estimates must equal the mean of the real values.

The second criterion states that the mean of the squared deviations must be minimal, which means that when the cloud of estimated values is more scattered compared to the cloud of real values, the estimate will be more inaccurate.

Methods

Depending on the stochastic properties of the random field and the various assumed degrees of stationarity, different methods of computing the weights can be derived, i.e. different types of kriging are applied. The classical methods include:

• Ordinary kriging assumes a constant unknown mean only in the search neighborhood of.
• Simple kriging assumes stationarity of the first moment over the entire domain with a known mean:, where is the known mean.
• Universal kriging assumes a general polynomial trend model, such as the linear trend model .
• IRFk kriging assumes to be an unknown polynomial in.
• Indicator kriging uses indicator functions instead of the process itself to estimate transition probabilities.
  • Multiple-indicator kriging is a variant of indicator kriging that works with a family of indicators. MIK originally showed great promise as a new method that would allow global mineral concentrations or grades to be estimated more accurately. However, these advantages are outweighed by other problems inherent to the practicality of modeling, due to the large block sizes originally used, as well as insufficient resolution at mining scale. In this case, conditional simulation quickly becomes the generally accepted replacement.
• Disjunctive kriging is a nonlinear generalization of kriging.
• Lognormal kriging interpolates positive data using logarithms.

Ordinary kriging

The unknown value is interpreted as a random variable located at , as well as the values of the neighboring samples . The estimator is likewise interpreted as a random variable located at , resulting from a linear combination of the variables.

To derive the kriging system, for the model assumptions the following error is admitted in the estimate: at is declared as:

The two quality criteria mentioned earlier can now be expressed as the mean and variance of the new random variable :

Unbiasedness:

Since the random function is stationary, , the following constraint is observed:

To guarantee unbiasedness of the model, the weights must sum to one.

Minimum variance:

Two estimators can both have , but the spread of their mean value determines the difference in quality between the estimates. To find the estimate with minimum variance, we need to minimize.

* for a detailed explanation, see the covariance matrix

* where the literals stand for .

Once the covariance model or variogram is defined, or , valid over the entire domain of analysis of , then we can write an expression for the estimation variance of any estimator in terms of the covariance between the samples and the covariances between the samples and the point to be estimated:

From this expression, we can draw several conclusions. The estimation variance:

• cannot be quantified for any linear estimator unless stationarity of the mean and of the spatial covariances or variograms is assumed.
• grows when the covariance between the samples and the point to be estimated decreases. This means that when the samples are farther from, the estimate becomes worse.
• grows with the a priori variance of the variable . When the variable has lower variance, the variance is smaller at any point of the domain.
• does not depend on the sample values. This means that the same spatial configuration (with the same geometric relationships between the samples and the point to be estimated) always reproduces the same estimation variance at any location within the domain. Thus, the variance does not measure the uncertainty of the estimate produced by the local variable.

System of equations

Solving this optimization problem (see Lagrange multipliers) leads to the kriging system:

the additional parameter is a Lagrange multiplier used in minimizing the kriging error to satisfy the unbiasedness condition.

Simple kriging

Simple kriging can be viewed as the mean and envelope of Brownian random walks passing through the data points.

Simple kriging is mathematically the simplest but the least general. It assumes that the expectation of the random field is known, and it relies on the covariance function. In most applications, however, neither the mathematical expectation nor the covariance is known in advance.

Practical assumptions for applying simple kriging:

• wide-sense stationarity of the field (stationarity of variance).
• Expectation everywhere zero: .
• Known covariance function

System of equations

The kriging weight in simple kriging has no unbiasedness constraint and is given by the simple kriging system of equations:

This is analogous to the linear regression of against the other .

Estimate

Interpolation by simple kriging is defined as follows:

The kriging error is given by the formula:

which leads to a generalized version of the Gauss–Markov theorem for least squares (Chiles & Delfiner 1999, p. 159):

Properties

(Cressie 1993, Chiles & Delfiner 1999, Wackernagel 1995)

• The kriging estimate is unbiased:
• The kriging estimate honors the actually observed value: (assuming no measurement error)
• The kriging estimate is the best linear unbiased estimator ofif the assumptions hold. However (e.g. Cressie 1993):
  • As with any other method: if the assumptions do not hold, kriging can turn out poorly.
  • There may be better nonlinear and/or biased methods.
  • When an incorrect variogram is used, no properties are guaranteed. However, a "good" interpolation is usually still achieved.
  • Best is not necessarily good: for example, in the absence of spatial dependence, the kriging interpolation is exactly as good as the arithmetic mean.
• Kriging provides as a measure of accuracy. However, this measure depends on the correctness of the variogram.

Applications

Although kriging was originally developed for applications in geostatistics, it is a general method of statistical interpolation that can be applied in any discipline to sample data from random fields satisfying the appropriate mathematical assumptions. It can be used wherever spatially related data (in 2-D or 3-D) have been collected and estimates of "gap-filling" data are needed at locations (spatial gaps) between the actual measurements.

Today, kriging is used in a variety of disciplines, including the following:

• Ecology
• Hydrogeology
• Mining
• Natural resources ^{[10] [11]}
• Remote sensing ^[12]
• Real estate valuation ^{[13] [14]}
• Integrated circuit analysis and optimization ^[15]
• Microwave device modeling ^[16]
• Astronomy ^{[17] [18] [19]}

Design and analysis of computer experiments

Another very important and rapidly growing area of application in engineering is the interpolation of data obtained as response variables of deterministic computer simulations,^[20] for example finite element method (FEM) simulations. In this case, kriging is used as a metamodeling tool, i.e. a black-box model built from a designed set of computer experiments. In many practical engineering problems, such as formwork design, a single finite element simulation run can take several hours or even several days. It is therefore more efficient to design and run a limited number of computer simulations, and then use a kriging interpolator to quickly predict the response at any other design point. For this reason, kriging is very often used as a so-called surrogate model, implemented within optimization procedures.^[21]

See also

• Bayesian linear statistics
• Gaussian process
• Interpolation
• Multivariate interpolation
• Nonparametric regression
• Radial basis function
• Regression kriging
• Space mapping
• Spatial dependence
• Variogram
• Gradient-enhanced kriging (GEK)
• Surrogate model
• Information field theory
• Lagrange multipliers